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Paper Discussions 


BE of the most important reasons for pre- 
senting a paper at a technical convention is 
that interested members of the audience 
can discuss the subject of the paper with the au- 
thor. Although these discussions may benefit the 
entire audience and add considerably to the value 
of the paper, it is unfortunate that discussions do 
not always evolve because the paper may not be 
well understood in a short presentation time or be- 
cause the interested members of the audience may 
not have time to formulate their question sen- 
sibly. Consequently, with a lack of discussion, 
many well organized sessions are not as valuable 
as they should be. 

This problem has been alleviated somewhat at 
the WESCON Convention and other technical 
meetings on the West Coast. Here fewer papers 
are presented in each session, but a panel is se- 
lected to discuss the papers after thev have been 
presented. Audience participation is encouraged, 
and the panel may include the authors themsel ves 
in a round-table discussion after all the papers 
have been given. Since the panel is given some 
previous knowledge of the subjects to be presented, 
a discussion does occur, and, subsequently, mem- 
bers of the audience contribute their thoughts 
which add considerably to the value of the paper 
presentations. 

There are other methods of obtaining discus- 
sions of papers, and one of them is to obtain writ- 
ten discussions which are published in an appropri- 
ate journal. This is a particularly significant 
method of obtaining discussion of TRANSACTIONS 
ON AUTOMATIC CONTROL papers because, although 


the papers receive the benefit of reviewer’s com- 
ments, they are not always presented at a confer- 
ence, and they are not usually exposed to open dis- 
cussion. We are fortunate to have a few discus- 
sions of papers in this issue. We believe that the 
most valuable discussions are made voluntarily, 
and we would like to receive more discussions for 
publication if they add to the clarity or broaden 
the scope of any paper either directly or by com- 
parison with an alternate method. 

The only other requirement, which may lead to 
some editorial censorship, is that the discussions 
by both the discussor and the authors must be 
written in a courteous manner. Of course, it is not 
always easy for an author to reply politely to a dis- 
cussion which appears to be inappropriately or un- 
justly made. However, this could be a fault in the 
organization and clarity of the paper itself, not 
obvious to the author, which could have misled 
many readers, and any further clarification, how- 
ever unnecessary it may seem to the author, could 
add to the value and prestige of the paper. 

We hope that the few discussions in this issue 
will be found interesting, and that they will evoke 
more discussions in future issues. Our correspond- 
ence section will be open for notes, sometimes brief 
papers, and discussions, as well as announcements 
and all informal documents not subject to review. 
Actually the correspondence is published as rap- 
idly as possible without the review delay, because 
we believe that more of this informal correspond- 
ence coupled with the larger, formal papers will 
add considerably to the interest and value of the 
TRANSACTIONS. —The Editor 


263 


te 


264 


IRE TRANSACTIONS ON AUTOMATIC CONTROL 


September 


The Issue in Brief 


Stability of Systems with Randomly Time-Varying Parameters, A. R. 
Bergen 


Due to intermittent failures, a system may alternate in a random 
fashion between different configurations having stable and unstable 
modes. The resulting stability of such a system is discussed in this 
paper. It is investigated conveniently by using matrix notation and 
Kronecker products. This paper is of particular interest because it is 
an excellent example of the powerful mathematical techniques that 
are becoming more widely used in the analysis and synthesis of control 
systems. 


Dynamic Programming Approach to a Final-Value Control System 
with a Random Variable Having an Unknown Distribution Func- 
tion, M. Aoki 


A number of papers have been written in recent years about the 
application of dynamic programming to control systems. In this paper 
the technique is applied to a final-value control system. Given a cri- 
terion of performance, the deviation of the process, and the domain 
of the control variable, a sequence of control variables is determined 
as a function of the state vector of the system and time in a manner 
that optimizes system performance. 

It is shown that the determination of the control variables must 
be made sequentially rather than specified as a function of the initial 
state and time because of the stochastic nature of the problem. By 
means of a functional equation technique of dynamic programming, 
a recurrence relation of the criterion function of the process is derived. 

The concept of suboptimal policy is introduced and numerical 
experiments are performed to test suboptimal policy suggested from 
a numerical solution of the recurrence relationship. 


Mathematical Aspects of the Synthesis of Linear Minimum Re- 
sponse-Time Controllers, E. B. Lee 


The synthesis of minimum response time control systems has 
been of interest for many years. In this paper a procedure is presented 
for finding, systematically, the forcing function to be applied to the 
process to give time-optimal control. The results are limited to the 
case where the process to be controlled has dynamics which are com- 
pletely known and are adequately described by a system of linear 
differential equations. 

The implication of the method presented is that any process as 
limited above can be time-optimally controlled provided that a 
system of transcendental equations can be solved. The results have 
been used to indicate that the method can be applied in quite general 
situations. 


Statistical Design of Digital Control Systems, J. T. Tou 


Although the problem of statistical design of digital contro] sys- 
tems has been considered in the literature, a simplified design pro- 
cedure is introduced in this paper. The useful minimum mean-square 
error performance criterion is adopted as the performance criterion, 
and the system design is based upon the Wiener-Kolmogoroff theory 
of optimum filtering and prediction. The treatment of this optimiza- 
tion problem makes use of the modified z-transform technique, and it 
is presented in close parallel with the procedure for the statistical 
design of linear continuous-data control systems. By making use of 
some important properties of the auto-correlation sequence, the cross- 
correlation sequence, the pulse spectral density and the pulse cross- 
spectral density, it is shown that the statistical design of digital con- 
trol systems can be carried out in a simple and systematic manner. 


Determination of Periodic Modes in Relay Servomechanisms Em- 
ploying Sampled Data, H. C. Torng and W. E. Meserve 


Sampled-data feedback systems containing a relay and a zero- 
order hold have been extensively investigated. The describing func- 
tion and the phase-plane approach has been used with some success 
in studying these systems, but these methods have been rather un- 
wieldy. 


In this paper, the difference equation approach is used to reveal 
all possible periodic modes of relay sampled-data feedback systems of 
any order in a simplifed way. This is achieved by expanding periodic 
functions of an integral argument in terms of orthogonal functions 
in the discrete domain. 

The effect of the presence of a dead zone in the relay characteristic 
is studied. The steady-state response of the system corresponding to 
step, ramp, or sinusoidal inputs is also discussed. 


Analysis of Pulse Duration Sampled-Data Systems with Linear 
Elements, R. E. Andeen 


Most of the extensive literature concerning sampled-data systems 
is based on a pulse amplitude modulation process. In this paper, the 
sampled data is represented by constant amplitude pulses, with the 
pulse width modulated by the data. The characteristics of these sys- 
tems are described, and an analytical means for studying the per- 
formance of these systems is presented. Analytical results are com- 
pared with experimental tests, and a short discussion of the relative 
advantages and disadvantages of pulse duration modulation in con- 
trol system is included. 


The Analysis of Cross-Coupling Effects_on the Stability of Two- 
Dimensional, Orthogonal, Feedback Control Systems, D. B. Newman 


Multipole or multidimensional control servos have been receiving 
considerable attention in the past few years, and this paper should 
serve as an excellent introduction to the subject for those not ac- 
quainted with the problems involved and the methods of analysis. 
In addition to a development of an algebraic model of the system, 
matrix methods of analysis are introduced, and the performance of 
the system is interpreted in terms of the pole-zero plots of the system. 
A technique showing how an unstable system may be made stable by 
the introduction of appropriate cross-coupling is discussed. 


Stabilization of Linear Multivariable Feedback Control Systems, 
E. V. Bohn 


This paper is another in the field of multipole control systems. 
Like the previous paper, it may serve as an introduction to the sub- 
ject although it is somewhat more general. This paper is also con- 
cerned with stabilization techniques for multivariable feedback con- 
trol systems. A method of system stabilization is discussed which re- 
duces the problem to essentially that of a single variable system for 
which stabilization techniques are well developed. This method uses 
a compensating matrix which is either the inverse system—or trans- 
posed system-matrix. An example of this method of system stabiliza- 
tion is given and its usefulness discussed. 


Correspondence 


A paper by G. J. Murphy and N. T. Bold! has resulted in an inter+ 
esting discussion between J. B. Cruz and one of the authors, Murphy. 
A note on a particular nonlinear feedback system is provided by B. 
Chatterjee, and Z. Bonenn presents an illuminating comparison be- 
tween various methods of determining stability including the describ- 
ing function and Lyapunov’s second method. The root locus may be 
modified to factor mth order polynomials as described by J. D. Glomb, 
and methods of determining the transfer function coefficients of a 


linear dynamic system from frequency response characteristics are 
discussed by P. V. Rao and E. L. Bai. 


PGAC Membership Directory 


This list of control engineers is again published to denote the 
active members in the field. Ze 


_ 1G, J, Murphy and N. T. Bold, “Optimization based on a square-error criterion 
with an arbitrary weighting function,” IRE Trans. oN ANTENNAS AND PROPAGA- 
TION, vol, AC-5, pp, 24-30; Jan., 1960, 
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Stability of Systems with Randomly 
Time-Varying Parameters’ 


A. R. BERGENT, MEMBER, IRE 


Summary—tThe stability of a system which alternates between 
stable and unstable configurations at random times may be investi- 
gated conveniently using Kronecker products. The system is stable 
with probability one if, and only if, all the eigenvalues of a specified 
matrix lie within the unit circle. If stability in the presence of a param- 
eter adjustment is to be investigated, a root-locus plot of the eigen- 
values is convenient. 


I. INTRODUCTION 
\ CONSIDERATION of systems with randomly 


time-varying parameters by Rosenbloom! in 

1954, and more recently in two papers by Ber- 
tram and Sarachik? and Samuels’ indicates a growing in- 
terest in such systems. This paper is concerned with 
random systems of a particular type, those in which the 
system alternates between two known configurations at 
random times. The alternation -between configurations 
is a purely deterministic process; the system remains in 
each configuration for a random interval. It is assumed 
that the random intervals are mutually independent 
random variables. 

One may presume such systems to arise as a result of 
an intermittent faulting of the system in some way, by 
an intermittent component failure or by a random inter- 
ruption of a control or guidance loop. In either case it is 
assumed that the system adopts one of two possible con- 
figurations. For simplicity it is assumed also that the 
order of the system is the same in either configuration. 

In a case of interest to control system designers one 
of these two configurations is unstable. The question 
then arises as to the stability of the system when this 
unstable state is intermittent. Presumably the system 
will be stable if the unstable state is induced infre- 
quently and is of short enough duration. 

Since the system switches randomly between stable 
and unstable modes it is necessary to define precisely 
what is meant by stability. A time-invariant linear sys- 
tem, initially in equilibrium, is said to be stable (strictly 
stable) if, and only if, the system returns to equilibrium 
after any finite disturbance. In the case of the random 

4 


* Received by the PGAC, January 13, 1960; revised manuscript 
received, May 26, 1960. This research was supported by the Elec- 
tronics Res. Directorate of the AF Cambridge Res. Center, Air Res. 
and Dev. Command, under contract AF 19(604)-5460. 

+ University of California, Berkeley, Gait, 

1A. Rosenbloom, “Analysis of Linear Systems with Randomly 
Time-Varying Parameters,” Proc. Symp. on Information Networks, 
Polytechnic Inst. of Brooklyn, Brooklyn, N. Y., vol. 3, pp. 144-154; 
April, 1954. ; De She : 

2]. E. Bertram and P. E. Sarachik, “Stability of circuits with 
randomly time-varying parameters,” IRE TRANS. ON Circuit THE- 
ory, vol. CT-6, pp. 260-270; May, 1959. ~ ; 
~ 3 J. C. Samuels, “On the mean square stability of random linear 
systems,” IRE TRANS. ON Circuit THEORY, vol. CT-6, pp. 248-259; 
May, 1959. 


system presently under discussion, however, the be- 
havior following a disturbance can only be described 
probabilistically. 

The return of the system to equilibrium is analogous 
to the convergence of a sequence of random variables, 
and, borrowing the terms of the statistician, the system 
may exhibit convergence in probability, convergence in 
the mean (norm), or convergence with probability one. 
Of these, convergence with probability one is the most 
desirable type of behavior assuring that for any bounded 
disturbance the system returns to equilibrium after suf- 
ficient time with probability one; this behavior matches 
closely that of stable time-invariant linear systems. For 
the purposes of this paper a stable system will be de- 
fined as one which returns to equilibrium with proba- 
bility one after a bounded disturbance. 


II. ANALYSIS: 


The problem under consideration is the stability of 
the autonomous randomly time-varying linear system 
described in the &th interval by the vector differential 
equation, 


ey'= Aw), baste 


x(t) is an n-vector called the state of the system with 
components x(t), x2(t), «+ + , Xn(t), called the state vari- 
ables; for an mth order system there are 1 such state 
variables. In the kth interval, A; is a constant nXn 
matrix. 

The solution of (1) is* 


a(t) = eAn(t—te-v y(t, 1), bay th 


HE NI qty (2) 


this being the solution in a particular (the kth) interval. 
In cases of practical interest the matrices A; are 
bounded, and therefore the state variables are contin- 
uous functions of time. Then 


x(th) ee eAkTky(t,_1) k= is Dr Cpe, (3) 


where 7),=te—t-1 is the kth interval. 

Eq. (3) relates the system state at the end of the kth 
interval to that at the end of the k—1st interval and, by 
iteration, to the initial state x(to). The stability of the 
system depends on the behavior of x(t) after a large 
number of random transformations, the typical one 
given by (3). It should be recalled that the randomness 
of the system resides in the intervals 7); the A; are a 


4B. Friedman, “Principles and Techniques of Applied Mathe- 
matics,” John Wiley and Sons, Inc., New York, N. Y., pp. 122-124; 
1956. 
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purely deterministic sequence of alternatively stable and 
unstable matrices.° 

Eq. (3) is convenient for investigating the limiting 
mean of «(t;), and this will be done to illustrate the tech- 
nique to be used later when stability is considered. 


Iterating (3), 


a(ty) = eAbTkeAb-iTh-1 . . . GARTRY(f5) | (4) 
and taking expectations (the intervals are assumed inde- 
pendent), 


E{ «(t,)} = Ej Atte} By eAbiTe-s} te Ej 4171} x(t). (5) 


x(to) may be assumed to be an arbitrary nonrandom 
disturbance. Assuming now that the alternate intervals 
are identically distributed random variables, alternate 
terms of (5) will be identical, and for Rk even, 


E{ x(t,)} = [E{ 64272} Ef eit} |k/24(t5) 


Wal?x(to), (6) 


where V4 is an abbreviation for the bracketed expres- 
sion. For k odd (6) is premultiplied by Ej ein} ; since 
this single transformation cannot affect the limiting be- 
havior of x(f,), it is sufficient to consider (6) with k even. 

From (6), the méan of x(t,) tends to zero with if, 
and only if, the eigenvalues of Wy lie within the unit 
circle. This condition is analogous to the condition im- 
posed upon the poles of a stable sampled-data system, 
1.e., the poles must lie within the unit circle. 

The problem of evaluating Y4 may be simplified con- 
siderably if the eigenvalues of 4; and A, are distinct. In 
this case, if the integral converges,® 


E{ eAxte} -{ et Tp. (TdT, = Mi(— Ax) 
o pee Uae 7) 


Here p,(7;) is the probability density function of stable 
(unstable) intervals. M,(—A;,) is a matrix obtained in 
the following way:’ find M;(s), the Laplace transform 
of p.(T%), and replace s by — A; and 1 by the identity 
matrix. Note that it is not necessary for the distribution 
of intervals to be the same in the stable and unstable 
cases. 

The convergence of the mean of x(t.) with & does not 
imply stability since it is entirely possible for the mean 
to converge while the variance diverges. It is therefore 
necessary to investigate the behavior of a norm of 
x(t). It will be shown, however, that it is possible to in- 
vestigate the behavior of a norm in a manner entirely 
similar to that used for the mean. This may be accom- 


> Clearly the stability of the system does not depend on whether 
A, is stable or unstable. Assume then that it is stable. 

° To test the convergence of the integral replace the matrix A; 
by the real part of the eigenvalue with the largest real part. If the 
integral does not converge, the mean of x does not converge:and it is 
therefore unnecessary to compute Wy. In the case of the matrix 
A, with eigenvalues having negative or zero real parts, the integral 
always converges; the test is unnecessary. 

7 Friedman, op. cit., pp. 113-117. 
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plished by the use of Kronecker or direct products of 
matrices. 

Kronecker products are related to tensor products; 
for a treatment of the subject the reader is referred to 
MacDuffee,* or Bellman.® Bellman!" has used the Kron- 
ecker product in problems related to that treated here. 
Kronecker products of matrices are defined in Appendix 
I, An important property of such products needed in the 
derivation which follows is also given in Appendix I. 

Consider the behavior of the vector y(t), which is the 
Kronecker product of the state vector x(t) with itself. 


y(t) = x(t) @ «(8). (8) 


The cross symbol denotes the Kronecker product. Con- 
sider y(t) in a particular interval: 


i) = #(1) ® x(t) + x) ® #(0) 
A,x(t) @ x(t) + x(t) @ Azx(t) 
[(4, @1 + 2 ®@ Ai) ly. (9) 


The last line in (9) follows by premultiplying x(t) (in the 
second line) by the identity matrix, using (15) of Ap- 
pendix I, and (8). Under consideration then is a new 
system (with state vector y) which has been derived 
from the old system (with state vector x). The mean of 
the new vector y may now be investigated in the same 
way in which the mean of the x vector was investigated. 
But since the vector y includes among its components 
the squares of the components of the x vector, it is pos- 
sible to investigate the mean-square behavior of the 
state variables. 

The use of (9) without modification leads to unneces- 
sary computational labor. The vector y contains among 
its elements components such as x1%2 and also xox1; these 
elements are identical and it is possible, by eliminating 
such redundant elements, to reduce the order of the 
matrix equation from dimension zn? to dimension 
n(n+1)/2. How this reduction proceeds may be illus- 
trated as follows. Suppose the same element appears in 
the second and fifth rows of the y vector. Eliminate this 
element the second time it appears, 7.e., eliminate the 
fifth row. At the same time modify the matrix which 
premultiplies y by adding the fifth column to the sec- 
ond and then eliminating the fifth row and column. Con-. 
tinue until all the redundant elements are removed. 
Clearly, this process removes all the redundant elements 
while maintaining the correct relations between the re- 
maining elements. 

Let the reduced matrix and vector be designated by 
B;, and 2 respectively. Then (9) after reduction is 


a(t) = Byz(t), thor St < te k=1,2,--- 


I 


(10) 


§ C. C. MacDuffee, “The Theory of Matrices,” Chelsea Publish- 
ing Co., New York, N. Y., pp. 81-86; 1956. 

* R. Bellman, “Introduction to Matrix Analysis,” McGraw-Hill 
Book Co., Inc., New York, N. Y.; 1960, 

* R. Bellman, “Limit theorems for non-commutative operators, 
I,” Duke Math. J., vol. 21, pp. 491-500; September, 1954. 

UR. Bellman, “Kronecker products and the second method of 
Lyapunov,” Math. Nachr., vol. 20, pp. 17-19; May, 1959. 
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with solution (at the switching times) 


@(ty) = ePePkz(t, 1), Em bee) esac (11) 


Eqs. (11) and (3) are now in the same form. The 
method previously used to investigate the convergence 
of the mean of x(&) can now be used to investigate the 
convergence of the mean of 2(t,); in the derivation 
simply replace x by z and A by B. It follows that if the 
eigenvalues of B,; and By, are distinct” and the integral 
corresponding to (7) converges, the behavior of the 
mean of 2(f,) depends on the locations of the eigenval- 
ues of 


Ve Ma(— B.)M,(— Bi): (12) 
If these eigenvalues lie within the unit circle, the mean 
of z tends exponentially to zero and so do each of its 
components which include E(x’), E(x”), +--+, E(xn?). 
In this case the expected norm for the disturbed system 
tends to zero with k. Moreover, because of the expo- 
nential convergence it is possible to establish the 
stronger condition of convergence with probability one. 
This is shown in Appendix III. 


Example 


The stability of the feedback system shown in Fig. 1 
is to be determined. The minor loop shown opens inter- 
mittently; the system is stable with the loop closed and 
unstable with the loop open. Assume that the loop is 
closed during the first interval. 


r(t)=O + 
oe, 
wy 


Fig. 1—Random system of example. 


The scalar differential equation describing the autono- 
mous system in the Ath interval is 


¢+(-1)F+¢=0. 


Choose for the state variables, x1; =c and x2=¢. Then in 
the &th interval 


Ae ae 
te) LA -(—1)*4 aed’ 
and in the notation of (1), 


2 In Appendix II, it is shown that while the eigenvalues of the 
_matrix involving the Kronecker products [shown in (9)] are not dis- 
tinct, those of the reduced matrix By ordinarily are if the original 
matrix A; has distinct eigenvalues. In any case the eigenvalues of the 
reduced matrix B, may be enumerated using the rule of Appendix II. 


Fae te 
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2607 
oe | 0) 1 | 
cal as Nie 
Consider next the matrix 
0) 0) 1 0) | 
ane ores  ) 0) 1 | 
0M 1) eee 
Ona 0) ie | 
0 1 0) 0 
Sea ty Oe a0 
0 0 0) 1 
0 0 —1 oie 
0 1 1 0 
—1 (-—1)* 0 1 
—1 0 (—1)*} 1 
Q  —1 —1 . 2(—1)*} 


which is found by using (13) in Appendix I. The y vec- 
tor corresponding to x is 


y= \ a1”, 1X2, Lox, aca?} , 


Here the second and third rows are seen to be equal. 
The third row,may be eliminated, yielding the z vector 


B= { x12, X12, ac9?} 


with the corresponding B; matrix 


hela, 0 
B= ai (pe 1 
f) Dua argh Ge a eae 


The B, matrix is formed from the preceding matrix of 

Kronecker products by adding the third column to the 

second and then eliminating the third row and column. 
In particular then 


0 2 0 
W Sie 1 lle | 1 1 
(ae 2 2 
0 4, 0 
Bos} —l +1 1 
0 22 


Consider now the eigenvalues of B; and By respec- 
tively. These can be computed directly or, as shown in 
Appendix II, they can be specified in terms of the eigen- 
values of the lower order matrices Ai and A» respec- 
tively. The eigenvalues of A; are 


those of B; are therefore —1+jV/3, —1, £1 9r/3. 
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Similarly, the eigenvalues of A, are 


those of B, are therefore 1+j+/3, 1, 1—7/3. Since the 
eigenvalues of B; and By are distinct, (12) may be used 
to investigate stability. 

First, however, a model for the random opening of 
the minor loop must be chosen. 

Assume, for example, as a plausible model, a purely 
random opening and closing of the minor feedback loop. 
Let the average interval during which the loop is closed 
(stable) be fixed at 1 second while the average interval 
during which the loop is open (unstable) is adjustable; 
let this average interval be 7». 

In either case, the probability density function of in- 
tervals is exponential 


1 ty 
DAT) = Sj e Tht 


with Laplace transform 
M,,(s) = (Ts + il Feat 


T, is the average interval. 

If Ts is sufficiently small the system presumably is 
stable, while if TZ, is large enough, system instability is 
to be expected. System stability as a function of 7, will 
now be investigated. 

Note first that the defining integral for /.(— Bz») does 
not converge for JT, greater than unity. Then for Ty, 
greater than unity the system is unstable. For T. less 
than unity the eigenvalues of Vz as given by (12) must 
be investigated. For the present model 


Vez =o (—T,B, + I)(— B, + vee}. 


Equivalently, and more conveniently, the eigenvalues of 
Wz! may be examined. In Fig. 2 are shown the root-loci 
of the eigenvalues of Vz as T. is varied from zero to one. 
For T, less than approximately 4, the system is stable. 
For T»2 greater than approximately 4, one of the eigen- 
values lies outside of the unit circle and the system is 
unstable. 


CONCLUSION 


In general it is difficult to obtain strong conditions 
defining the stability of systems with randomly time- 
varying parameters. In the particular case considered 
here it is possible to obtain sharp stability conditions, 
1.e., the system is stable with probability one if, and 
only if, all the eigenvalues of a specified matrix lie 
within the unit circle. The system for which this applies 
is one which alternates between stable and unstable con- 
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Fig. 2—Root-loci of eigenvalues. 


figurations. The times of switching from one state to the 
other are random. 


APPENDIX | 


The Kronecker product of two nXn matrices A and 
B (denoted by A ®B) is an n? Xn? matrix defined as fol- 


lows: 
° ee 
Ann B ; 
ai; are the elements of A. 
In particular the Kronecker product of the state vec- 
tor x with itself is the 2? dimensional vector 


an1B ay.B-:- 


13 
GniB \ 


4@B=| 


: F Poet 
xOx= { x12, Bia, ee hs Bae oka, eae ees 


X2Xn,° *° ) XnX1, XnX1, pe an}, (14) 


which includes among its components the squares of 
the components of x. 
An important property of Kronecker products is 


(AB) @ (CD) = (A @C)(B ® D). (15) 


APPENDIX II 


Under consideration are the eigenvalues of the re- 
duced matrix B,. Assume that the eigenvalues of the 
original » Xn matrix A; are distinct. Then 


Ax&i = wk: 


where the &; are the m linearly independent eigenvectors 
corresponding to the distinct eigenvalues u;. Consider 
the right hand Kronecker product of (16) by Jé;. Using 
(15), 


(16) 


(A; @ T(E ® &) = wilE: @ &). (17) 
Consider also the left hand Kronecker product of 
Ape = mjk; (18) 


by Jé; 
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(I ® Ax) (Ei ® &) = uwj(E; @ &). (19) 
Adding (17) and (19), 
[4, @I+1®@ Al: @&) = (atu) @é). (20) 


The matrix on the left of (20) is just the matrix in (9), 

e., the B, matrix before contraction. The eigenvectors 
of this matrix may be identified as the n? Kronecker 
products of the m eigenvectors of the matrix A; having 
for eigenvalues the n? possible sums of eigenvalues of the 
original A; matrix. Included among these sums are 
terms such as #i+pe and we+pi; hence the eigenvalues 
are not distinct. 

Consider now the eigenvalues of the reduced matrix 
B,. First note from (20) that 


[4, @7+1® A, |[(E @ &) + (& ® &)] 
= (ui + w)[(E: @ &) + (E&; @ ED]. (21) 


Now £;@£;+£&;@&; is a vector which displays the same 
structure regarding redundancies as does x@x. This 
can be shown by writing this vector as a sum of self 
Kronecker products. 

This being the case, a contraction of the type dis- 
cussed in the text of this paper eliminates n(m—1)/2 of 
the redundant elements leaving the relation between 
the remaining elements intact; 7.e., the vectors which 
are the reduced versions of £;@£;+£;@&; are the eigen- 
vectors and the w;+y; are the corresponding eigenvalues 
of the matrix B;. But now, it is to be noted, the inter- 
change of indices does not yield a new eigenvector; the 
n(n—1)/2 eigenvectors which were generated by such 
interchanges have been eliminated. The eigenvalues of 
B,, therefore, are all possible u;+p;, with 7 less than or 
equal to 7. 

While the eigenvalues of A; are distinct, it does not 
follow that those of B, are distinct, although this will 
usually be the case. For example, if the eigenvalues of 
A; are —1, —2, and —3, those of B, will be —2, —3, 
—4 —4, —5, — 
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APPENDIX III 


Assume the eigenvalues of Wg lie within the unit 
circle. Then the infinite sum of vectors E{z(t,) } con- 
verges: 


oO 


wat Ej 2(ts)} < %. 


(22) 


The necessary and sufficient condition for this con- 
vergence is that V,* tend to zero with k;" this is assured 
if the eigenvalues of VW, lie within the unit circle. Con- 
vergence of the sum is also assured for each component 
of E} (tk) \_ In particular 


a 1x2(th)} < © EOE LAG, wee 23) 


But now, by Tchebycheff’s inequality, 


Ms 


P{| x(t.) | > 2(h)} < ©, 


= Dale 
ke} k=1 


where € is an arbitrary positive number. Then by the 
Borel-Cantelli lemma x,(f) can exceed the bound e 
only a finite number of times in an infinite run. Since 
this holds for arbitrarily small e, x,(t,) tends to zero with 
probability one. Since this is true of all the components 
of x(t,), x(t,) itself tends to zero with probability one; 
the system is stable. 
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Dynamic Programming Approach to a Final-Value 
Control System with a Random Variable Having 
an Unknown Distribution Function” 


MASANAO AOK 


Summary—As described in the introduction of this paper, the 
state vector x,, representing the present state of a sampling control 
system, is assumed to satisfy a difference equation 


XLn+1 = Ti Gtr, ln; Un) 


where v,, is the control vector of the system subject to random dis- 
turbances r,. The random variables r,, are assumed to be independ- 
ent of each other and to be defined in a parameter space in the fol- 
lowing manner: Nature is assumed to be in one of a finite number, 
q, of possible states and each state has its own parameter 
value, thus specifying uniquely and unequivocally the distribution 
function of 1,. It is further assumed that we are given the a priori 
probability z=(zi, - - - , 2.) of each possible state of nature, 


q 
ete aE etd ys eee ineay & 0: 


t=1 


Given a criterion of performance, the duration of the process and 
the domain of the control variable, a sequence of control variables 
{y,} is to be determined as a function of the state vector of the 
system and time, so as to optimize the performance. 

The sequential nature of the determination of v, is evident here, 
since the stochastic nature of the problem prevents specification of 
such a sequence of v’s as a function of the initial state and time. By 
means ot the functional equation technique of dynamic programming, 
a recurrence relation of the criterion function of the process, k,,(x, z), 
is derived, where x is the state variable of the system when there 
remain n control stages. 

When q is taken to be two and the criterion of performance to be 
xy? with the constraint on the control variables 


N=" 
> #2 SK, 


i=0 


where xy is the final state vector of the system, it is shown that 
k,(x, z) is the form w,(z)x? and the optimal v,,(x, z) is linear in x, as 
might be expected. The dependence of k,(x, z) and v,(x, z) on z is in- 
vestigated further, and w,(z) is found to be concave in z. Explicit 
quadratic forms for w,(z) are obtained in the neighborhood of z=0 
and 1. The optimal v,(x, z) is found to be monotonically decreasing 
hi 

When the domain of v, is restricted to a finite set of values, as in 
contactor servo systems, no explicit expressions for k,(x, z) and 
v,(x, z) are available. However, k,(x, z) is still concave in z and ap- 
proximately given by 


zhn(w, 1) + (1 — 2)Rn(x, 0). 


By solving the recurrence relation numerically, this approxima- 
tion is found to be very good for moderately large n, say 10. This 
means that if one has explicit solutions for =p, and p=», then 
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one can approximate the criterion function for the adaptive case 
where Pr(p= p:) =z by a linear function in z as shown above. 

This provides fairly good lower bound on &,(x, z), and in turn 
serves to determine an initial approximate policy. 

A concept of suboptimal policy is introduced and numerical ex- 
periments are performed to test suboptimal policy suggested from 
numerical solution. The system behavior under the suboptimal pol- 
icy is also discussed. 


I. INTRODUCTION 


A. Mathematical Descriptions of General 
Control Systems {2] 


ET us denote by Sa dynamic system of » degrees of 
ile freedom. Then, a set of scalar functions of time 
i, x(t), ¢=Ip2,2--, 2 can describe completely 
the state of S at each time instant, 7.e., given x‘(é), 
i=1, 2,---,m for t<to, the behavior of S with time 
for {>t is completely described. These functions of 
time can be regarded as m components of a vector func- 
tion of time x(t), called the state vector of the system S. 
Although in a most general case the value of x(t) for 
t>t) is dependent not only on the current state vector 
x(to) but also on all the past values of the state vector 
x(t), £<to [1], in the following only such systems whose 
future state vectors x(t) are completely and uniquely 
determined by their current state vector x(to) =c are 
considered. That is to say, under suitable existence and 
uniqueness conditions, x(t) becomes a function of time 
t and of the initial condition vector c, x(t)=<x(c, 2), 
x(0) =c, taking tp =0 without loss of generality. 
From the uniqueness property which x(t) satisfies, 
there follows 


aC AL Se to) oa x(x(c, ti), to), (1) 


which is nothing more than a mathematical representa- 
tion of the principle of causality [2]. 
Let us take a unit time interval At and set 


ol Cope AN eka Oo (2) 


1.e., the initial state vector c is transformed into a new 
state vector 7(c), after the lapse of a unit time interval. 
Then from (1), 


wl, MAL) = TLC * Tc) too) Pe) a) 


which is to say, the state vector at t=nAt is given by 
T"(c), the mth iterate of the function 7(c). This means 
that the time behavior of the deterministic dynamic 
system S at time t=nAt is determined by the successive 
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iterates of a specific function of x, or successive point 
transformations on c [2]. Here, the kind of transforma- 
tion that is applied on the state vector at t=nAt, to get 
a new state vector at the next time instant f=(n+1)At, 
is independent of x(mAt). On the other hand, feedback 
control systems are designed in such a way as to utilize 
more than one type of transformation as a function of 
the state vectors [3], [4]. 

Let ti, t2, - - - be a sequence of times, #;<f2< --- at 
which the system S is subject to a change in the form of 
a transformation, then 


«(tep1) = T(x(te), r), (4) 


where the control vector v; is included to indicate the 
fact that the transformation is dependent upon the 
choice of the control vector. 

The choice is made in such a way as to optimize some 
performance index assigned to the control system 
[5], [6]. 

For example the state vector of a linear sampling con- 
trol system in which the control vector enters linearly is 
given by 


C1 = Ad Bie Pri; (5) 


where 


At is the sampling time interval, 

x; =x(RAt) is an J-dimensional state vector, 
A is I by I matrix, ; 

Bis I by J matrix, 

v, is J-dimensional control vector at t=kAt, 
r, is [-dimensional vector at t=kAt. 


The term ~% may represent the effect of a random 
forcing variable with Markovian property [7] on the 
system. Then the relation of (5) denotes a stochastic 
transformation on the state vector xz. ; 

More generally, the state vector at time (w+1)At is 


given by 
0, Ue 2, See, (6) 


the type of transformation on x, being governed by vari- 
ables v, and 7,.! 


Xn+1 = Ete, Uny Ti), i= 


B. Adaptive Control Systems 


In many control systems, there exists random effects 
of one kind or another. They may be due to noisy com- 
ponents in the systems or due to a noisy environment 


which disturbs the outputs of the system. Very often, 


we do not have complete information about the distri- 
bution functions of random variables. Thus in many 
control processes, decision-making devices or controllers 
are called upon to make decisions under various uncer- 
tainties. As the process unfolds, additional information 


1 Eq. (6) means that once the system is in the state «p, regardless 


~ of how the system reaches the state xp», the next state xn41 1s given 


by (6). The way of looking at the next state and an optimal control 
vector as functions of current state vector, and not as functions of 


- the initial vector, is essential in dealing with nondeterministic control 
systems {4], [6]. 


”] 


as 
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may become available to the controller which may shed 
light on these uncertainties. The control system then 
has the possibility of “learning” to improve its per- 
formance based upon experience, 7.e., the decision-mak- 
ing device may “adapt” itself to its environment. 

In this paper, a class of distribution functions is con- 
sidered in connection with certain time-discrete adap- 
tive control systems, using a mathematical model of 
Bellman and Kalaba [2]. Distribution functions for 
random variables are assumed to be members of a cer- 
tain parameter space, which is assumed to be discrete 
and finite. Therefore, a random variable R has the dis- 
tribution function of the form F(r, a), where 


PCa = a;) = 


q 
S oo ite 
41 


as (Oy. Pema yD eras ng 


Let us define the states of nature to be determined by 
these g possible values which the parameter @ can as- 
sume. The state of nature is assumed not to change with 
time. Therefore, a control system must be able to esti- 
mate the state of nature in which it is likely to be 
operating. 

When 2;=1 for some 7, then the control system is 
known to be in the 7th state of nature. That is to say, 
in this case, it is known that the random variable R has 
the distribution function F(r, a;). This situation will be 
referred to as stochastic to distinguish it from the fol- 
lowing situation. If 0<z;<1 for all7=1, 2, - - - , g, then 
the control system is called adaptive. 

The fact that this situation is really adaptive in the 
sense described at the beginning of the section can be 
seen as follows: 

In the adaptive situation, there is a positive proba- 
bility that the random variable R has the distribution 
function F(r, a), 7=1, 2,--+-,4¢. 

Since an optimal choice of control vector v at each 
time instant depends on the knowledge of the. distribu- 
tion function of R, here the system is called upon to 
produce control vectors without having exact knowledge 
on the distribution function of R. It is also assumed 
that the a priort probability z; is modified after each 
control action, as will be shown later, to produce a pos- 
tertort probability 2,’,7=1, 2,---,q. That is to-say, 
the current estimate on the parameter @ is updated at 
the conclusion of each control action. 


II]. FINAL-VALUE SYSTEMS: GENERAL DISCUSSION .. 


Let us take as the performance index of a control sys- 
tem some function ¢ of the state vector xy at the time 
instant ty, at which time the control action is assumed 
to terminate. The time ty is called the final time. Control 
systems with such a performance index are called final- 
value systems [8|-[10]. The control system S must 
decide to apply v which minimizes the estimated ex- 
pected value of (xy), given the present state variable x, 
the number of the remaining decision stages ” and the 
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present estimate of the probability that the system is 
operating in the ith state of nature, 2,,7=1,---, @. 

The problem is adaptive, because the control device 
must guess or learn the state of nature and act accord- 
ingly, and this estimate as to the actual state of nature 
must be revised at the conclusion of each decision stage 
depending on the history of the observations of the ran- 
dom force r. 

For the ease of presentation, let us take the Bernoulli 
random variable 


(are with probability p; 
— \=c with probability 1 — p;,7 = 1, 2 
Pr (p = pi) = 2, Pr (p = po) = 9 — 2, pi > Pe. 


That is, p is now the parameter a of the distribution 
function F(r, x) and nature is assumed to be in one of 
two possible states, H; and H». The random force r is 
assumed to be independently and identically distrib- 
uted for each of the WN stages. 

Let us define the following: 


(7) 


Tn 


h,(x; p:) =the expected value of ¢(xy) employ- (8) 
ing the optimal policy, given the pres- 
ent state variable x of the system, 
and the number of the remaining 
decision stages ”, and given that the 
state of nature is the 7th state, z=1, 2. 


If we know definitely that nature is in the 7th state 
(i.e., z=1 or 0), then the optimal control force at each 
stage is given by the v, which is obtained by solving the 
recurrence relation (9) derived by the application of the 
principle of optimality [2]. 


I(x; p:) = min [po(at) + (1 — p)o()] 
h,(“; pi) = min [piltnr(at; ps) + 1 — Pi) Mn—1(%-; pi) | (9) 


n-1 


t=, 350 > 
jae 


ee NG 


where xt and x~ are the state variable of S at the next 
decision stage, when the present state variable is x, with 
ry=-+c and r=—c respectively. 

When we are given only the a priori probability as to 
nature being in the state H; or He, (9) must be rewrit- 
ten by defining ,(x, 2) to be the quantity corresponding 
to h,(x; p:) given the present estimate z that the nature 
is in H;. Noting that the a posteriori probability at the 
nth stage becomes the a priori probability for the next, 
(n—1)st, stage, 
ki(x,2) = min {z[pid(at) + (1 — pio(o-)] 


+ (1 —2)[polat) + (1 — po ]}, 
fin(w, 2) = min {2[prkaa(at, 2!) + 1 = pr)hrales2”)] 


On-1 


+ (1 — 2)[pokn—a(at, 2’) 
Bike (1 wy p2)Rn-a (xe, 2'’)]}, (10) 
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where 2’ and z” are given by (11) and (12) respectively. 

If the present estimate of nature’s being in H, 1s 2, 
and if r=c is realized, then the a posteriori probability 
that nature is in H, will become 


z 1 
ae ‘ \pige aed de 
vA = 2) P92 ame 
pi p (eee 
Z 
When r= —c is observed at this stage, the a posteriori 
probability becomes 
2(1 — pi) 1 
eRe) meaner To 
2(1 — — z)(1 — po — 2% 
. : 1-f- ao 
z 


where a and a, are the likelihood ratios of H,; over He, 


after r,=—c and 7,=-+c are observed, nameiy 
l= po pe 
ayo = ) a=—-: 
LPP Pi 


After » such observations, the @ priort probability 2 
becomes 


(13) 


An 
z 


where a, is the likelihood ratio of the ” observations, 
1.e., if r=c are observed 7 times, and r= —c are ob- 
served 7» times, then 


\%/1 — po\™ 
o = (2) ( *) ; N+ ne =n. 
1! a Reogeey 


From (10), 


(14) 


Ri(x, 2) = min { [261 + (1 — 2) pol d(x) 


+ [s(t — pr) + (1 — 2)(1 — pr) |o(a)} 
ka(x,2).= min.{ [zp + (1.— 2) pelenur(xt,.2’) 


On-1 


+'[2(1 = p1) + (1 42)(L — p2) |e G29} 
is B03 ee PS 


Let us note that if z=1 or z=0; then 2’=2’’=1 or 
z’=32''=0, and (10) and (15) reduce to (9), since 

If one could obtain, before solving (15), some informa- 
tion on the functional structure of k,(x, z) from the 
knowledge of h,(x; p1) and h,(x; p2), then one would be 
in a better position to devise an appropriate computa- 
tional procedure and/or an analytic approximation of 
solving (15). 

One such structural information is given by the fol- 
lowing. 
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Proposition 
R,(x, 8) of (10) and (15) is concave in 2, i.e., 
Rn(x, 2) > 2k,(x, 1) + (1 — 2)Rn(x, 0), 
m= 1,2 


aS (16) 
For the proof, see Appendix I. 

Eq. (16) supplies the lower limit on &,(x, 2), given 
Rn(x, 1) =h, (x, Pa) and k(x, 0) =hn(x, p2). 

The upper limit on &,(x; 2) is given from (10) by 


R(x, 2) < min max [pika_a(xt, 2’) + (1 — pdRnala-, 2”), 


Un-1 


Prine", 2) + (1 — ps) Raila, 2’) 


min max [k,_-1(x*, 2’), kn—1(a-, 2”) | 


Pn-1 


IA 


Clyla, &). 


Observe that ct,(x, 2) requires less computation than 
fenke,) 2). 

Thus, we know that the minimum of the estimated 
expected value of $(xy) for this type of adaptive final 
value system, k,(x, z), is concave in gz, and (16) provides 
a lower bound on k,(x, 2), given the expected values of 
the criterion function for the corresponding stochastic 
case, h,(x; p1) and h,,(x; po). 

Let us note that in order to prove the concavity (16), 
the actual form of the difference equation of the system 
(6), distribution function F(r, a) and ¢ are immaterial. 
The concavity is a general characteristic of final value 
control systems witha finite number of states of nature. 
Since the knowledge of a sequence of optimal v is equiva- 
lent to that of the criterion functional values [2], the 
lower bound of &,(x, z) may be used to derive initial 
approximate policy to start a sequence of approxima- 
tions in policy space [11]. 

Let us assume for the moment that the recurrence 
relation (10) has been solved, and let us consider cer- 
tain fixed values of x and n. Then in (10) let us define 


Ss; and Ss as 
ee pln ae, 2 (1 Piao 2 pt 1, 2,--(18) 


thus s; and s, will be in general functions of x, n, v and z. 
We can now write (10) as 


min [zs + (1 — 2z)se| 


I 


kn(x, 2) 


= 25;* + (1 — 2) 50%, (19) 


where s,* and s9* are s; and sy for a certain v which is 
optimal for x and considered. 

If one regards s; and sz as the loss of the control sys- 
tem with v when nature is H, and HZ respectively,” 
then the control process may be considered to be the S 
game [12] where nature has two strategies. When na- 
ture employs a mixed strategy with probability distri- 
bution z=(z, 1—2), the expected value of the loss to the 


2 This is exactly true only for z=1 or z=0. 


va 
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control system with v becomes 
2s, + (1 — 2)s50, 


and the s* =(s1*, so*) is the minimum of (20). 

Viewed in the light of game theory, the use of the a 
priort and the a@ posteriori probability distributions in 
(10) seems to be quite natural. This can be formalized by 
assuming that z will be transformed by the Bayes 
formula. 

The reason for introducing the concept of S game 
here is that there exist mathematical theories on set- 
theoretical relations among the classes of strategies, and 
they are useful in discussing optimal strategies or opti- 
mal policies. Although it does not seem possible to fit 
questions in adaptive control processes completely in 
the existing frame of the theory of statistical decisions 
and sequential analysis, certain analogies can be used to 
advantage in the construction of the theory of adaptive 
control processes, and in deriving approximate solutions 
for functional equations derived by the application of 
the principle of optimality. 

When g, 2’, and 2” are close to each other, it is reason- 
able to expect that k,(%,-2),/k,(xy 2) and k(x) 2") are 
also close to one another, since they are continuous in 
x and zg. The recurrence relation defining &,(x, z) by 


hi(x, 2) = kala, 2) 
k(x, 2) = min {2[Pikn—(x*, z) + (1 — piRn—a(x-, 2) | 


+ (1 — 2) [pokn—i(at, 2) + (1 — po) Rn_s(a-, 2)]} 
itd 3: ane 


may be used as an approximation to the exact recur- 
rence relation (10). 


(20) 


III. StrucTURE DEPENDENCE OF QUADRATIC CRI- 
TERION FUNCTION OF FINAL VALUE CONTROL 
SYSTEMS ON THE PARAMETER 2. 


Let us take a quadratic form for ¢(«wn) as 
o(xw) = (aw’, Px), 


where P is a positive semi-definite matrix, xy is the final 
position of the system, and «y’ is the transpose of xy. 

Let us assume that the system is to operate under the 
constraint 


(22) 


N-1 
Di (mn’, Don) < K, (23) 
n=0 
where D is a positive definite matrix and v, is the control 
variable of the system at the mth decision stage and »,,’ 
its transpose, and it is assumed that —» <v,<+0, 
n=0,1,--:, N—1. 

The problem can then be formulated to minimize the 
expected value of 

N-1 
In = (n', Paw) + 27 (mn, Din), d > 0, 


n=0 


where d is the Lagrange multiplier [6], [13], [18]. 


(24) 


ro 
~I 
HS 


The criterion function of the system is given by 


min 
D0sOls sere oN 


fy (Xo, 2) = HIN) (25) 


where xo is the initial position of the system, z is the a 
priori distribution of the states of nature between A, 
and Hs; and vp, v1, - - - , ¥y_1 are to be chosen optimally, 
z.e., to realize the minimum value of A(Jy) taken with 
respect to the random force r and g. 

Let us assume further that if the state variable of the 
system at the present decision stage is x, then the state 
variable at the next decision stage will be given either by 


te= Ax a bg, orby 24. = Ax + BB», 6) 


where suffixes + and — refer to the two possible states 

of the random variable r given by (7) to which the sys- 

tem is subjected. Here, the effect of 7 is, therefore, as- 

sumed to appear in two possible B’s, 7.e., By and B_. 
For ease of presentation, let us assume here that the 

state variable of the system is one-dimensional, and 

carry out the discussion in terms of scalar quantities. 
Eq. (24) is rewritten as 


N—-1 
In! = In/p =n? +d Do mn’, (24’) 
n=0 
where \d/p is renamed as \, \>0. 
The criterion function is now taken to be 
fy(%o,2) = min E(VJy’), (25’) 
00; VL 50.05% oN 1 


since the minimization of (24) is equivalent to that of 
(24’). 

The new state variable of the system is related to the 
old state variabies by 


t, = ax-+ b4v, and x_ = ax+ b_2, (26’) 


where 0<a<i. 
The functional equation is given similarly to (15) by 


fx(x,2) = min {[zp. + (1 — 2)polfra(xs, 2’) + [A — pi)z 


Speer lie) par v=,.2") 1 Noy a? 
he 1, ig es 
i (422). min { [zp1 + (1 — 2) po] (aa + b400)2 + rv0? 


vo 


+ [(1 — pile + (1 — pe)(1 — 2)](ax + b_m)2}, (28) 


where 


eng pis eee (1 — pi)z 
piz + (1 — 2)p2 (1 — pie + (1 — pe)(1 — 2) 


From (28), (0f1/0 v0) =0 gives the minimizing vo. The 
existence of the minimum is assured by the fact that 
f(x, 2) is convex quadratic in vp. 

This vo is given by 


(27) 


Yoo go(2) “x, (29) 
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where 
go(z) = — aE(b)/(A + E(O’)), 
E(b) = 2£\(b) + (1 — 2) F2(b), and 


E,() = pbs + 1 — pb, 


and E(b?) is defined similarly to E(6). 
The criterion function is given by 


i= 1,2. 


(30) 


fil Uy. 3). = Wl sexe 
where 


wi(z) = a*[1 — E(b)?/(A + E(b?))]. (31) 


It can be shown inductively that, in general, v, is 
linear in x and f,(x, z) is quadratic in x and can be ex- 
pressed by 

fr(x, 2) = Wn(z)x*, 
Val, 2) = Sn Zee (32) 


The recurrence relations for w,(z) and gn(z) are given 


by 


E(bwn(z))? | 
A + E(b?w,(z)) 
= AND ea) Vee le 


SRO a*| Blws(2) = 


(33) 
E(bw,(2)) 
yh Bb) 


gn(z) = — ES Py ei ee 


The symbols appearing in (33) and (34) are defined as 
follows: 
E(wn(z)) in 2E1(Wn(z)) ree 2) F2(Wn(z)), (35) 


where 


Ej(wn(z)) a PiWn(2’) + (1 ha: PjOne’), i= i} 2 


Similarly, 
E(bw,(z)) = 2Ei(bwn(z)) + (1 — 2) E2(bwa(z)), (36) 

where 
Ex(bwn(z)) = pibywnlz’) + (1 — pi)b-wa(z”), i = 1, 2; 
and 

E(b’wn(z)) = 2Ei(b°w,(z)) + (1 — 2) E2(b?wn(z)), (37) 
where 
E(b'wWn(Z)) = pidy?wn(2’) + (1 — pi)b?w,(2”), i = 1, 2. 


Now that the recurrence relation on w, is given, let us 
investigate the functional structures of g,(z) and w,(sz), 
namely their dependence on g. Assume 


O< poSpi<1,b_ <b, and B2<d,2, (38) 


the case where b_>6,, can be easily taken care of in the 
similar manner. 


1960 


Eq. (29), when written out in detail, becomes 


go(z) = — a| Eo(b) + [E1(b) — Bo(b)]2}/{x + £2(b2) 


+ [£1(6?) — E2(b?) Jz}, (29’) 
where 
Ex(b) — E2(b) = (pi — p2)-(b — 6) > 0. (39) 
Similarly 
E,(b?) — E2(b?) = (pi — po) (by? — 6?) > 0. (40) 


Since 


(=) 

dz 

Exe) = Ea(0) | [A+ 2(b*)]— E2(0) [Ex(0*) ~ Ex(0*) | 
{A+ Ea(b?) + [Ex (62) — Eo(62) |}? 


does not change sign for 0<z<1, go(z) is monotonic in z. 
For J sufficiently large,’ go(z) is monotonically decreas- 
ing. 

From (31), 


wi(z) = a?| 1 — [£2(b) + (Ei(6) — Eo(b))z]?/[\ + E2(b2) 
+ (£,(b") — E2(b?))z]}. (41) 
If one specifies b4 to be 


b4. = Sm (42) 


then, since E;(b?) =s?,7=1, 2, wi(z) becomes a concave 
quadratic function in gz. 

Or, even without assuming (42), if X is large, then 
w,(z) will be concave quadratic in zg, neglecting terms 
0(1/d?).4 

In the following, let us consider the situation where 
E(b?)/A<1 holds.* Loosely speaking, this means that 


3 It is true that with equality in (23), the value for will be 
uniquely determined. 

In the dynamic programming approach to the use of Lagrange 
multipliers, however, it is more convenient to work with (24) or 
(24’) directly, treating \ as if it is another independent variable, 
rather than working with K in (23). (This is particularly true in 
computational solutions.) 

For example, f, and v, are tentatively determined by assigning a 
certain value to X. Inserting these v, back in (23), if it is found that 
(23) is not satisfied, then the value of » is increased and the above 
process will be repeated. Thus \ can usually be interpreted as price 
or cost. 

Parameter studies on K can be replaced by equivalent parameter 
studies on \. The assumption on the range of K can be replaced by 
an equivalent assumption on the range of d. This is what 1s done in 
the paper. : : 

The theoretical justification for the above statement is required 
and can be found in Ch. 4 of [18]. 

4 The symbol 0(x) means a quantity of the order x, thus 0(1/d*) 
indicates a quantity which approaches zero essentially at the rate 
of 1/X?, 4.e., 4 


— constant asA— ~. 


6 Asymptotic expansions in terms of 1/) are therefore investi- 
gated. - 


pa 


Aoki: Dynamic Programming Approach to a Final-Value Control System 


no 
~I 
Or 


only a small amount of energy is available for control 
purposes, since \ can be interpreted as price of control 
as mentioned in footnote 3. 


From (41), 


1 
wi(z) = aft — < | Eo(b) + (£,(b) — F())3F 


+ 0(1/22). (43) 

As a first approximation, take 
Wits) Ge 1, (44) 

Then, from (33), 

a (e)a== aa (45) 

This, in turn, gives from (34) 

qontl 

gu(2) = — —— | Ea(b) + [Zi() — Ex(0)Jz}. (46) 


Therefore, g,(z), and consequently the optimal v, are 
monotonically decreasing in g. Note that g,(z) =0 for 
the zero-th degree term in the expansion in 1/)X. 

Let us make (45) more precise by retaining terms of 
the order up to 1/d. 

Then, 


1 
wi(z) = oft Pat [£2(b) + [E£i(b) — r(6)|5' ese 
and 
1 
wnra(e) = | Beals) — —- BGw(3))*], 
t= 2, or N = 1.2 445") 
As can be seen from (35) and (37), the right-hand 


side of (45’) contains w,(z’) and w,(z’’). They are ex- 
panded in Taylor series at 2z,° 


eve @ +() @ - 
*Wr\S = Wr\S a ee VA 


1 /0?wy 
+2(S8) @ t+ 0@ — 99, 


a) Oe" ) 
Go Le 
02 /, 


Z 


Wn(2') =< Wn (2) a= ( 


il Cae, 
+>( 2) (el! — 9 + O((0" = 9), (47) 


Then, from (35) 


Fine) = ee (=) Z, (48) 


02? 


6 At z=1, the derivative from the left is used. At z=0, the deriva- 
tive from the right is used. 
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where 

2Z = (pi — p2)?22(1 — 2)?/{ [zp1 + (1 — 2)p2] 
LEE Spa Gh eG ai, 


Z & (f1 — p2)?2"/2p2(1 — pe), if zis close to zero. (49) 
If z is close to one, then 
Zp (Py ape) (leis) 27 2 Pi Vi) (490) 


Following the similar derivation, 
g 


OWy, 


Oz 


B(bun(2)) wn) 20) + ( ) a1 = )[2206) — F206) 


O2Wn 
5 Ze, 
02? 
where 


22’ = [pz + (1 — 2) po|b4(2” — 2)? 
+ [(1 — pire + (1 — po)(1 — 2) ]b_(2” — 2)? 


(50) 


When 2, 2’ and 2” are close to each other, Z and Z’ 
can be approximately made zero. 
From (43’), (45’), (48) and (50), 


1 
We(z) = aft nae [E2(b) + (E1(b) — E2(6))z]? 


+ 2(Ei(6) — E2(6))-Z + otB(0)*} 


Since in the neighborhood of z=Oand z= 1, Z isquadra- 
tic in 2, W2(z) is quadratic in z in the same neighborhood. 
Also w,(z) will be quadratic in z, in the same neighbor- 
hood as can be seen inductively, if the terms up to 1/ 
are retained. 

Differentiating (45’) twice and making use of (48) 


and (49) 
0? O7Wn, 
a(t oa) Ga) 
02? 02? 


~ 2a!(E,(b wee 
moe! ( 1) — Bay)? (51) 


O07? Wnt 1(2) 
O22 


1.€. 


O7Wn41 
Oz? 


07Wp 
— a(t + Ao.) ( ) 
0,1 02" Jo,1 
Jqint2 


r 


— 


(E1(b) — £2(b))?, 


pup)? 


= TOG ghee 
0 po(1 — po) 


pil = pi) 
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Since 
0? 2a? 
ee eee 
02? aN 
0° We 2a‘ 
== (E,(b) — E2(b))?[(1 + Ao) + 2]. 
02? 0,1 mn 


In general, 
2a 
r 
[po,® + p91" -t 4. 2 a] 


= — ee (0) — BW)? 


OP Wr4 


02? 


0,1 


(Ei(b) — E2(d))? 


(52) 


where 
Pt 1 = Ao,1. 
Therefore, in the neighborhood of z=0, 1, 


(a?po,1)” 


(E1(b) — E2(6))? 
a? \7tl 
ae 
Po,1 
3 s 
(e) 
P0,1 
+ Wn41(0)(1 — 2) + Wnyi(1)z. 
If z=1, then zg’ =z’’=1. Therefore, from (35)—(37), 


Wnti(Z) = 


2(1 — 2) 
(53) 


E(wn(1)) = Er(wa(1)) = prwn(1) + (1 — pr)wn(1) 
=a (1). 
E(bwr(1)) = Ei(bwn(1)) = (p1by + (1 — pi)b_)wi(1) 
= E1(b) wp (1). 
E(bwn(1)) = Ey(b?wn(1)) = (pbs? + (1 — p1)b_*)w,,(1) 
= F,(6")w,(1). 
From (33), 


1 
Wn41(1) ~ a bed Ex(0)tan(t) | 061) 


n=1,::-,QN, 
wi(1) = a*[1 — E,(b)?/(A + E,(6?))]. 


Similarly, if z=0, then 2’=z’"’=0 and 


(54) 


Wn+1(0) ~ ae = | 00), 
n=1,---,QN, 


wi(0) = [1 — B2(6)2/(A + Ea (82))]. (55) 
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From (32) and (53), 
Sn(z) ~ — a/d-[E(b)wn(1)z + E(b)wn(0)(1 — 2) 
+ 2(1 — 2)(Ei(6) — E2(6))- (w(t) — wa(O))] 
+ 0(1/22). (S6) 
Since its derivative with respect to z does not change 
sign as seen from 
dgn 
dz 


= — (a/d)[(Ex(d) — E2(b)) wn (0) 


+ E1(6)(wn(t) — wn(0))], (57) 


£n(Z) is monotonic in zg. 
More generally, if terms of only up to 1/X are re- 
tained, since 


E(bw,(z)) = E(b)a* + O(1/A), (58) 
Gn(z) = — a2*1/)[Fo(b) + (E£1(b) — Eo(b) E1(b) E2(6))z] 
+ 0(1/A?), Gy) 


g,(z) is monotonically decreasing in z. 

From (32), the optimal v,(x, z) is, therefore, mono- 
tonically decreasing in zg under the conditions stated 
above. 

Thus, f(x, 2) is quadratic in the neighborhood of 
z=1 and 0 and concave over 0<z<1. 


IV. ContTActTor SERVO FINAL VALUE SYSTEMS 


Let us next consider the situation where a constraint 
on v is imposed as 


SIL, 


and let us consider the difference equation 


Lat = On + tn + Uy (60) 
Lae i with probability p (7) 
: —c¢ with probability 1 — 9, 
where p is assumed to be known. 
The criterion function is taken to be 
o(av) = xy’. (61) 


The state variable x, is taken to be a scalar. 


The functional equation describing the control proc- . 


ess is given similarly to (9) by replacing p: by ~, 
n(x) = min [pxy? + (1 — pa] 


v=tm 


(9’) 
(x) = min [phya(xy) + (1 — p)ina(x-)] 
o=t 
Ee 
where 
x, =ax+tce+v, x_=ax—c+2, (62) 


- from (60). 
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Since the equation governing the system is fairly 
simple, a certain amount of analytic manipulation of 
the criterion function is possible. 

From (9’) and (62) one gets 


poe + (ax + (2p — 1)¢ + 24)?].(63) 


Ai(x) = min [4p(1 


vu=tm 


Therefore, the optimal control force is given by 


Cee or 4-25 —.1)c = 0, 
De hae ee (64) 
and 
a) = Ape peaeet Op = 1yetooil4, (aso 
Similarly, 
ho(x) = [phi(wy) + (1 — p)hi(w_)]os (65) 
where 


ty = ax 6 -— dp; t= CH = 6 3P wy, 


and v, is understood to take m or —m, whichever is the 
optimal choice in that it minimizes f(x). 
From (63’), 
hi(x4) = 4p(1 — pie? 
+ [ata + ac + ave + (2p — 1)¢ + 9,t]?, 
hi(x_) = 4p(1 — p)e? 

+ [a2x — ac + avot+ (26 — 1)e + o-]?, (66) 
where v;* and v; represent respectively the optimal 
choice of v1, when the initial position is x; and x_ re- 
spectively. 

From (65) and (66), 
ho(x) = 4p(1 — p)c?(1 + a?) 
+ [atx + avs + (2p — 1)c(1 + a)]? 
+ [o.? + 2pnit — (1 — p)orfac 
+ 2(a?x + ave + (2p == 1)c)(poit + (1 ao; por). (67) 
If vy+=v;-, this becomes 
ho(x) = 4p(1 — p)c?(1 + a?) 


+ fate + (2p — 1)c(1 +a) + ave+]%. (68) 


This will be the case if | ac| is sufficiently large so that 
the outcome of ry does not affect the optimal choice of 
v1. In other words, | «| is sufficiently away from the 
switching boundary for 71, so that two possible values 
of 72 can not make the system straddle the switching 
boundary. 

If vx = —v;+, then (67) becomes 


ho(x) = 4p(1 — p)c?(1 + a?) 
+ 8p(1 — p)acost + 4p(1 — p)(oi*)2 
+ fo2x + (2p — 1)c(1 + a) + ave + (2p — 1)0x+}?. (67/) 
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The value of v2 is determined by 


— msgn (a2a + (2p — 1)c(1 + a) + 21); 


He 
if =; 
ve = — msgn (ax + (2p — 1)c(1 + @) + (2p — 1)o1*); 
if v7 = — vt. (69) 


The general structure of hy(x) is seen to be given by 
h(x) = 4p(1 — p)[ce? + (ac + ait)? +--- 
+ (a6 + ev-1 aN —?ny_at)?| 


pnp aoe Peta) (leieetegre ct a". 1) Ys 


SOY or ie net Oe Ret oN on |?, (70) 
where 
‘¢ if v4 = 2,", 
é = 
11k v= — 0,7, 
| Vi; iO = pear 
Te ; 
(2) — shone Oi vy = — 957, 


and v;+ and v;- are defined similarly to v;+ and v;-. The 
truth of this formula can be seen by the application of 
mathematical induction. 

The control force is given by 


ww = — msgn(a%x+ (26 —1)citat:-:+a%—}) 
Sey a yy). (71) 


It should be noted here that although (70) gives the 
functional form of hy(x), the values of ¢;, y:, 7=1, 
2,:--, N—1 and vy are not given explicitly, and that 
(70) does not serve to compute /y(x) as the function of 
x and WN. 

Let us now investigate a certain dependence of the 
policy and the criterion function on x and p. The de- 
pendence of the % and /,(x) on p is expressed as %(x; p) 
and hy,(x; p), t.e., v(x”; p) will be either m or —m depend- 
ing on which gives smaller /;,(x; p) at the point (x, p). 

Arguing inductively, one sees 


Iy.(%; p) = bn(—x;1— p), &k=1,2,---,Nholds. (72) 


The proof is given in the Appendix II. 

The relation (72) is significant because it means the 
reduction by half of the amount of computation neces- 
sary for solving the functional equation (9’). This rela- 
tion holds, even when v is extended from v= +m to 
more than two possible values, so long as the values v 
can assume are symmetric with respect to v=0, as in 
v—m, 0, or —™m. 

The cases of y= +m and v= +m, 0 are computed and 
h(x) are found to satisfy M(x; p)=m(—x; 1—)), 
| uo ie A 

Looking at (70), one notices that the explicit depend- 
ence of the criterion function h,(x) on x decreases as n 
becomes large, and that once (x) becomes constant, 
then f(x) will remain the same for all n>k. 
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The recurrence relation (9’) was solved compu- 
tationally’? with the boundary condition that h,(x; p72) 
=h,(D; pi) for «>D and hy(x; pt) =ha(—D; p.) for 
x<—D, i=1, 2. The essential steps in computational 
procedures are sketched in Appendix III [6]. 

Having computed h(x; p1) and h(x; p2), one can use 
them to provide a lower bound on the criterion function 
of the adaptive final value system R,,(x, 2) of (10) as dis- 
cussed in Section II, where z=Pr+(p=hi), 1—2=Pr 
(p= pr). 

Fig. 1 shows &,(x, 2) as a function of g for various 
and x. It is seen that for »>8, the lower bound as a 
linear combination of h,(x; p1) and h(x; pe) is a fairly 
good approximation to k,(«, 2). 


kn(x,0) 108 


0.0170 0.3706 06294 0.8536 0.3706 0.6294 0.9830 


m =9/128 
-0,25<X<0.25 


kn (x,t) 108 kn(x,6):108 


7 5000 
4000 
3000 
2000 


1000 


n=8 


0.9830 01464 0.8536 


0.0170 0.3706 0.6294 0.8536 0.0170 0.3706 06294 0.8536 


Fig. 1—The criterion function, kn(x, z), as the function of 
the a priori probability z. 
V. ONE-STAGE SUBOPTIMAL POLICY 


Fig. 2 shows an optimal control variable as a function 
of x and a for p=0.625, D=4, a=, c=75, m=9/128, 


n<N=12.8 Notice that the boundaries between 
Vn(x; p) = +m and v,(x; p) = —m are not simple straight 
lines. 


Given the total duration N of the control process, the 
optimal sequence of control vectors {vn(x; p)} are 
given by solving (9’) sequentially for m=1, 2,---, N. 


7 The computation was done on the IBM 709 at the Western 
Data Processing Center, University of California at Los Angeles. 

§ Nis taken to be 12. For N>12, the system stays pinned to either 
x=D or x= —D, because of the particular boundary condition of the 
ee until the number of remaining stages m becomes less than 12 
or 11. 
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THE STATE VARIABLE OF THE SYSTEM, X, 


Xq7 70.25 Xp 770.125 by Xm 20 Xq 70125 Xp 70.25 
| eee e-----—---, --- ~~ 
° 
oie oe ace SS 4 
z 
ac Be r eat” Gn ee, 
zy *bF----H a 
zg r 
g 5 2535S Se —-------—+ 
oh 6 ee — — ts i 
oe SS et Ee 
a3 alt 
Se G32 s5 Se Se 
At. nl| 3 oa ee ee 
3 frat 
y 10 RP----- a 
= : 
i 2 SS SS 
12 
Ixal £0; O=1/4 — — —(x,n) WHERE m IS THE OPTIMAL CONTROL FORCE 


d=-(2p-i)c/a (x,n) WHERE —m |S THE OPTIMAL CONTROL FORCE 


d=-1!/56 FOR a=7/8, C=I/I6, AND p= 5/8 m=9/128 


Fig. 2—Boundaries between +m control force and —m control 
force for the optimal policy for p=0.625. 


When N=1, the optimal control vector is given by 


v1(x; p) = — m-sgn (ax + (2p — 1)c). (64) 


If (64) is used as v,(x; p) for n> 2, then it constitutes 
a suboptimal policy. 

The switching boundary, then, is independent of 
and is represented by a straight line given by 


x= — (2p — 1)c/a (73) 
and shown in Fig. 2 as the line x, =d. That is, according 
to this one stage suboptimal policy, v,=—m if x,>d 
and v,=-+m if x,<d. 

From Fig. 2 it is noted that for x >0 the optimal and 
suboptimal policies agree fairly well, although the agree- 
ment is not particularly good for x <0. 

It might be expected, therefore, that if the system 
starts from the initial position x =D, then the subopti- 
mal policy is a fairly good approximation for the optimal 
policy. 

To see to what extent the conjecture is correct, the 
Monte Carlo method [16] was used to simulate the sys- 
tem behavior from initial positions x» = + D. Forty sim- 
ulated runs were made using random numbers to gen- 
erate r,=-+c and 7,=—c with appropriate p. With 
i sD) 


E(2n?)optimai = 0.00367, 


} for p = 0.625. 
CON \eatoptimal SS 0.00374, 


With Xo = —D, 


WEEN )optinal = ae 
Ey) caboptinial ae 0.00441, 


In spite of the relatively similar values in E(xy’), the 
suboptimal policy for x» =D appears to be better than 
that for x)= —D, as conjectured, by the fact that out of 
40 trials with x) =D, the optimal and suboptimal policies 
gave the same xy? in 21 trials, but in 40 similar trials 
with «x)= —D, the optimal and suboptimal policies did 
not give the same xy? in any case. 

Although these results are by no means conclusive, 
they tend to support the conjecture that the one-stage 
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suboptimal policy is a fairly good approximation to the 
optimal policy for the control system under considera- 
tion. 

The adoption of this suboptimal policy simplifies the 
control problem considerably. 

The approximate analysis of the control system with 
the one-stage policy will be presented next. 

The state variable under the suboptimal policy is 
given by 


Xn41 = A%n + tn + In, (60) 
where 
Un = — m-sgn (ax, + (2p — 1)c), 
oe with probability p 7’) 
iS 
—c with probability 1 — p 


ORG ain O Ka <~— ile 

Let us- note thati 17, >0) thenx,01< ve wandeit 
Xn<Xe, then xn41>%, where x,= —(2p—1)c/a, because 
(m—c)>0 by assumption. Namely, if the system devi- 
ation from the origin (the position of the equilibrium in 
the absence of the random disturbance 7) is positive, the 
deviation at the next stage will invariably be smaller 
than the present one. If the present deviation is less 
than x,., then the system will move towards the origin 
at the next stage. 

HOT NAN oO 


Vn, = — M, 


and 


nei > Le — 6 — m= — (m+ 2c). 


That is to say, once x, becomes greater than 
x,=min [—(m+2pc), —(2p—1)c/a], then x, can not 
become smaller than x_, for all k>7, 1.e., x, is the lower 
bound. Similarly, one sees that 


xy, = m+ 2c(1 — p) 


provides an upper bound. That is to say, if x, <x,, then 
x,<x, for all k>n. Simulated system behaviors are 
functions of the sequence of random variable r,. That is 
to say, different orderings of +c and —c with 7 in {ta} 
will, in general, result in different behavior of x, with n. 
Fig. 3 shows three such typical simulated system be- 
haviors with one-stage suboptimal policy. 

One noticeable feature of the system behavior is that 
x, approaches the origin for a few stages consecutively, 
and then x, jumps away from the origin; and, from 
there, it again approaches the origin for the next few 
stages. 

Let us now investigate the number of stages before a 
jump occurs. This serves to show how the jerkiness of 
the system depends on system parameters a, c, m and p. 

Suppose that x, = —(m+2pc), and xo=x1. Let k be 
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the first time that x,>, is realized. The largest k will 
be realized when 7, is equal to —c for all k stages, and 
the smallest k will be given when 7,=-+c¢ occurs con- 
secutively. 

In the former 


(m — e)(1 — at) 


x, = — a*(m + Qpc) + 
(a4 
2p —1 
Saeed (74) 
a 
In the latter 
Gl k) 
L, = eller yee Claas ce 
mene OX 
2p — 1 
Sole comele (75) 
(64 
From (74) 
Be 2h 1)6 = 
log| a UP | log| m+ 2pe | 
aS 1—a a l=—a@ - (76) 
a log [a| 
From (75) 
(2p-1 + 
log [n+ p | log [m+ ape OS] 
1-a a —a 
k> BCNE, 


log [a] 


With a=, p=0.625, c=zy¢ and m=9/128, k>7 in (76) 
and k>1 in (77). Therefore, k may be expected to range 
from 1 to 8; k=8 is realized with probability (1—/)8, 
and k=1 is realized with probability p. 

If x» =x,, then the similar argument gives 
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Xy2 0.17 


a*0.875 
czI/i6 X"-0.148 
0.20 m?9/|28 Xo0.25 


Xe 


-0.20 


Fig. 3—Typical behaviors of the state variable xn, with the 
numbers of past control actions 1, for p=0.625. 


VI. CONCLUSIONS 


It was shown that for a class of final value systems 
considered in this paper, k,,(x, z) is concave in the pa- 
rameter 2 which represents the a priori knowledge on 
the states of nature. 


> 2p—1 2e 
log| @ =< - SP "| — tog | m + 2e(1 = 6) += "| 


1—a a 


The maximum b is given with — sign, and the minimum 
k is given with + sign in (78). The range of & is from 1 
to 11. The maximum & occurs with probability p*, and 
the minimum & with probability (1—)). 

Thus, after 2 becomes such that x,<x, <x, is real- 
ized, then the bounds on the number of consecutive 
stages k where the system approaches the origin mono- 
tonically are obtained. Since Fig. 3 shows x, starting 
from x9=D, and not from x)=x, or x,, the computed 
range of k need not be a good estimate for the range of 
k in this figure. As a matter of fact, the range of & in 
Fig. 3 is expected to be narrower, and this is indeed the 
case as seen from Fig. 3. 


log [a] 


(78) 


The structural dependence of k(x, z), when (xy) is 
quadratic in xy, was investigated in detail under the in- 
tegral constraint on the control variable v, and it was 
shown that kp(x, 2) =w,(z)x? and that w,(z) is quadratic 
in the neighborhood of z=1 or 0. The optimal control 
variable was shown to be monotonic in g, and linear in x. 

When the control variable is of the binary form +m, 
no explicit dependence of k,,(x, z) and a,(x, z) on « is ob- 
tained. 

However, one-stage suboptimal policy was shown to 
be a fairly good approximation to the optimal policy 
and was used to establish lower and upper bounds on the 
system deviations. 
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APPENDIX | 
THE CONCAVITY OF k,(x, 2) IN 2 


From (10), Ri(x, 2) is linear in z (special case of con- 
cavity). Assume that k,(x, 2) is concave in g. Then, 


Rn(xt, 2’) > 2/k,(at, 1) + (1 — 2’)Ra(at, 0), 
BAe, 8) 2 2 Ra(x, 1) + (1 — 2’) Ra(x-, 0). 
Also from (15), 
Rn+1(x, 2) 
2 min [zp1 + (1 — 2) po] [zh (at, 1) + (1 — 2’)Rn (xt, 0) 
ste aL a) TU 2) (1 — -p:)] 
- |e" Raa, 1) + (1 — 2’), (a, 0)]. 
Since from (11) and (12), 
[ep1 + (1 — 2)ps|2’ = zp, 
[2(1 — p:) + (1 — 2)(1 — po) Jz” = 2(1 — pr), 
using the relation 


min [f(s) + g(o)] > min f(») + min g(2) 


(79) 


(80) 


where f(v) >0, g(v) >0, (79) becomes 
Rngi(x, 2) > min [spik, (xt, 1) + (1 — 2) pokn(at, 0) 


a 2(1 "ea pdRA, 1) i (1 ce z)(1 iv p»)Rnlx-, 0) |, 
> z min [pik,(xt, 1) + (1 — prdka(a-, 1)] 


a (1 oy Z) min [pokn (a*, 0) “fe (1 x p2)Rr(x-, 0)| 


= Bnyi(x, 1) + (1 — 2)Rn4i(2, 0), (81) 
a 
Rnvi(x, 2) = 2hnpi(x, 1) + (1 — 2)Rn41(x, 0). 
Therefore, kn+1(x, Z) is concave in z. This completes 
the mathematical induction on n. 


(82) 


APPENDIX II 
Tue Proor oF h,(x; p) =h,(—x; 1—p) 
The relation (82) is assumed for all k. 
(x; p) = — %(—4;1— 9p). (83) 
Then 
I(x; p) = bn(—«; 1 — p) (84) 
- is true for all &. This can be seen as follows: 
For k=1, from (63’), 
h(a; p) = 4p(1 — pyc? + [ax + (26 — te + s(x; 9)? 
Ay(—x;1— p) = 4p(1 — pe? ; 
+ [-ax + (211 — p) — the + 1x(—2;1 — p)P? 
= 4p(1 — po? + [ax + (2p — 1c + r1(2; p)]?. 


Therefore, (84) is true for k=1. 


Aoki: Dynamic Programming Approach to a Final-Value Control System 


281 
Assume that it is true for all 2 less than nm. Then for 
k=n from (9) 


Tin(205 b) = [plina(at; p) + (1 — p)bnr(%-;P) loner, (85) 
Ten te Loop) = |(1i— ip ena (( 2) 1p) 


Te a a A ae icici, (te) 
where 
xt = at + ¢-+ n(x; p), x = ax — c+ 0,(x; Dp), 
(—4#)* = a(—«) +¢+7(—2;1— 6) = — ww) 


ty —%; 1 =p) “Pane; p) = — Ge) 
(=a) = a(—a) — ¢ + m(—a1 — p) = — @. 
Now, (86) can be rewritten as 
h,( =x; 1 — p) 
SAC = p)iiai(—( =a)", p) + ply = (2) 3) enters 
= [0 — p)hna(0-; 6) + Plena(0t; d) loner) 
= h,,(x; p). (87) 


This completes the mathematical induction loop on k: 
The truth of (83) can be proved in the following man- 
ner. 

For k=1, from (64) 


ox(x; p) = — m-sgn (ax + (2p — 1)o) 
n(— a; 1 — p) = — m-sgn (—ax — (2p = 10) 
=a go ep) (88) 
Assume that (83) is true for all indexes up to k. Then 
Vet 1(X3 p) 
= — m-sgen (ay 4+ (26 — 1)hed +a+---+ a") 


“eet Sh eee cate 
Pee — 4,1 — p) 


= — m-sgn[—a*tly — (2p — 1)cA+at--- +a) 
“ha sae iikales ye; 
where ¥; is y: for the (—x, 1—p) pair. 


Since 7; is equal to v,(«; p) or equal to (2b—1)2;+ as 
defined in (70), 
‘ ane +b) => ilepp) or 
Lap = Atowt(— a5 1 = 9) 
and therefore, ¥:;=—vy:i, t=1, 2,---, k, since 
v,;+(—x; 1—p) = —v;*(x; p) is true by assumption for 
1=1,2,---,k.’ This completes the loop of mathemati- 
cal induction on the index k. 


(89) 


(90) 


APPENDIX III 
COMPUTATIONAL PROCEDURES 


Let us consider a set C consisting of a finite number 
of x values from the domain of x, | | <D. Lét us 


denote elements of C by ¢1, C2, °° *, Cm. 
9 In the induction, it is assumed that u(x; p) = —v,(—x;1—p) 
holds for n=1, 2,--+ +, k. Hence, vn*(x; p)=—v,+(—x; 1—p) also 


holds for 7#=1, 2,+--, & 
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The functional value of h(x) is evaluated at each of 
x=c;,i=1,2,---,m. This can be done because explicit 
value for x, and x_ are given by (62) for each of x=c; 
as functions of v. 

The functional value of ho(x) at x=c; is evaluated 
next. It is seen from (9’) that to compute /e(c;), it is 
necessary to have /1(c14) and hi(ci_), where 


i: = OC oats @apP Ws 


Since in general ci, do not belong to the set C, it is 
necessary to approximate h,(ciz) from a set of func- 
tional values /y(c;), +=1, 2,---m. 

The same process is repeated for x=c;,1=2,-+.-.,™m 
to obtain h2(x) at each of the element of C. It is now ob- 
vious how to continue to obtain h,(x), n=3,4,---. 
Generally speaking, the larger the value of m, the better 
the accuracy of computation and the longer it takes to 
compute the solutions. 
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Discussion 


P. M. De Russo (Rensselaer Polytechnic Inst., Troy, 
N. Y.). The techniques of this paper are valuable, but 
they generally suffer from the fact that the systems 
which result are extremely complicated, and the prac- 
ticality of their physical realization is therefore often 
doubtful. Dr. Aoki says, “The suboptimal policy sim- 
plifies the control problem considerably and is a fairly 
good approximation to the optimal policy for the con- 
trol system under consideration.” But is this not saying 
that, in many situations, the conventional techniques 
which have been utilized for some time are more desir- 
able when practical physical realization is included in 
the specifications? It appears that the value of the 
techniques of this paper are that the “best” system is 
designed and can be utilized as a standard for evaluat- 
ing nonoptimum systems and for investigating the ulti- 
mate limitations to system performance. Also, I would 
like Dr. Aoki’s comments on the convergence of the 
discrete approximation in the continuous case. 

Author’s Comments: Dr. DeRusso pointed out cor- 
rectly one of the major uses of the techniques discussed 
in this paper, that is to say, the solutions of functional 
equations of dynamic programming give optimal solu- 
tions which serve as standards for evaluating nonopti- 
mal solutions. 

There exist, however, at least two other major ad- 
vantages of the techniques discussed in the paper over 
the more conventional techniques. 

One is the fact that dynamic programming approach 
furnishes not only optimal solutions, but also structural 
information on optimal solutions, which in turn may be 
used to advantage in devising “good” suboptimal solu- 
tions. It is to be noted that ordinary numerical pro- 
cedures will not give structural information on optimal 
solutions. Furthermore, the fact that dynamic program- 
ming admits a new kind of approximation, approxima- 
tion in policy space, which is a more natural approxima- 
tion in many realistic situations, cannot be ignored. 

Another point is the fact that, in spite of many restric- 
tive assumptions, the formulation of adaptive processes 
via dynamic programming, as done by Bellman and 
Kalaba, is the only general mathematical formulation 
of stochastic and adaptive control processes. 

In my paper, computational solutions are considered 
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only for sampling control systems, and only analytic 
approximations are discussed for the continuous control 
system of Section III. 

However, when the process is continuous, and it is 
desired to solve functional equations of dynamic pro- 
gramming technique numerically, it is necessary to use 
difference equations instead of differential equations. 
This is always the problem we are faced with in trying 
to use a digital computer in obtaining solutions to prob- 
lems which do not admit analytic solutions. 

The question of convergence of solutions of difference 
equations to those of corresponding differential equa- 
tions naturally arises. There are several papers [11], 
[14], [15], [17] in which this question is investigated. 
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The results, however, are still confined to rather special- 
ized type of functions. 

Although the dynamic programming technique is a 
powerful tool for dealing with multi-stage decision proc- 
esses, of which control processes are one special kind, 
actual solutions of functional equations are not always 
routine matters, especially when the dimension of prob- 
lems exceeds three or four. Computational solutions of 
problems of high dimensions are at present limited by 
capabilities of available digital computers, and by vari- 
ous analytic difficulties. This situation will be remedied 
as soon as more powerful digital computers become 
available, and when more attention has been focused 
upon the problem. 


Mathematical Aspects of the Synthesis of Linear 
Minimum Reponse-Time Controllers’ 


E. BRUCE LEE} 


Summary—A procedure is presented for finding, in a systematic 
matter, the forcing function to be applied to a process to give time- 
optimal control. The results are limited to the case where the process 
to be controlled has dynamics which are completely known and are 
adequately described by a system of linear differential equations. 
It has also been assumed that the process is not subject to unknown 
disturbances. 

The implication of the method presented is that any process as 
limited above can be time-optimally controlled, in those cases where 
a time-optimal controller exists, provided an underlying system of 
transcendental equations can be solved. This system of transcen- 
dental equations can be quite easily solved for systems of order less 
than three, even with two forcing functions, and for various other 
special cases of higher order. The results of a study which indicate 
that the method can be applied in quite general situations are being 
presented elsewhere [1]. 


HistoricAL BACKGROUND 
\ ROUND 1950 a number of people became inter- 


ested in the control improvement that could be 

achieved by the use of nonlinear feedback. One 
of these nonlinearities was the relay, where it was estab- 
lished, on the basis of plausible arguments and experi- 
mentation, that the relay controller, if properly 
switched, would give a smaller response time than any 
other controller subject to the same amplitude limita- 
tion. The first results were very restrictive as to what 


* Revised manuscript received by the PGAC, April 25, 1960. The 
material of this paper was first presented as “Synthesis of Minimum- 
‘Response-Time Controllers” at the IRE Mid-American Electronics 
Conf., November 4, 1959. ; ; ; 

+ Minneapolis-Honeywell Regulator Co., Minneapolis, Minn. 
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could be controlled, but there were many results, and 
all were of great engineering significance as evidenced 
by the results achieved [2 |-[16]. 

One of the first restrictions that had been imposed was 
to restrict the results to processes which had real dis- 
tinct characteristic roots. This first restriction was in 
part removed by Bushaw [17], where he considers sec- 
ond order processes with complex characteristic roots. 
See also [12] and [27]. Another restriction imposed was 
that the differential equation describing the process not 
contain derivatives of the forcing function; 7.e., if the 
plant was described by means of a transfer function, the 
transfer function could not have any numerator dy- 
namics. This problem was discussed by Smith [19], 
where he introduced a new coordinate to handle the 
jump involved when the relay switches. Hung and 
Chang also considered this problem [18]. A little survey 
of the literature shows that this problem had been en- 
countered in the statistical design of control systems. 
The general transformation to eliminate numerator dy- 
namics can be found in [20]; this result can also be found 
in Stone, ef al.} 

One problem inherent in the design of time-optimal 
controllers, as well as in all control concepts that con- 
trol the state vector, is that all elements of the state 
vector be measurable. For processes of order two or 
three, it may be possible to generate the state vector by 


1 See [27], p. II-7. 
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successive differentiation of a measured variable. In 
practice, it is hard to obtain derivatives of higher than 
two or three. This problem was attacked by Harvey,? 
where he shows that most systems of high order are 
basically coupled second order systems and that, in 
most cases, there does exist a transformation to measur- 
able coordinates. 

One of the first extensions of the concept of relay 
control to systems of a quite general nature was made 
by Bellman, Glicksberg and Gross [23] or [24], where 
they give a proof of the previously accepted theorem 
which states that relay control is time-optimal. Their 
developments were restricted to processes with real dis- 
tinct characteristic roots and to processes where the 
control matrix is assumed to be non-singular. 

Bass [22] presented without proof several important 
theorems connected with the optimum relay controller. 
He was the first to indicate that the differential equa- 
tion adjoint could be used in determining the zeros of 
the forcing function. LaSalle [25], [26] has also pre- 
sented a number of important theorems related to time- 
optimal control. Work in Russia, primarily influenced 
by the fine leadership of Pontriagin, has established 
quite general results in the area of time-optimal control. 
Russian results contain many important general items 
connected with the existence and uniqueness of optimal 
controllers. Of great general significance is Pontriagin’s 
principle of the maximum. See [28]-[33] for a com- 
plete account of the Russian work. The problem of 
time-optimal control has also been studied by Desoer 
[34] using classical calculus-of-variation techniques 
and by Krasovskii [35 |-[38]. Also of a related nature 
is the paper by Rozonoer [39], which discusses the use 
of Pontriagin’s maximum principle and its equivalence 
in many cases to Bellman’s dynamic programming [24]. 

In all this work, the most important question of en- 
gineering interest is: when does one switch the relay 
sign? In most of the results it is concluded that the so- 
lution to the adjoint differential equation will provide 
the zeros of the relays. The results of this paper show 
that the basic structure of the switching equation is con- 
nected with a set of simultaneous transcendental equa- 
tions, and that to find the relay switchings one must in 
some manner solve these equations. 


STATEMENT OF THE PROBLEM 


The process control problem of concern here is to 
select an allowable control vector a(t) (allowable as de- 
fined below) to carry the process, as represented by the 
state vector #(¢) from an initial state £ to a prescribed 
state, &(¢), in minimum time. 

If a(t), the forcing vector, is restricted to a closed 
region, Q, where specifically | #;|<m., i=1,2,---,r 
then the control vector is said to be allowable. 


’ 


* [bid., ch. 2. 
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The process under control to be considered here is 
assumed to be governed by the linear vector differential 
equation 


dx 
— = Aa) + BAL, (1) 
which has the solution [40] 
a(i) = &() bo + if BN EUN)a(d) a, (2) 


where £ is the initial condition vector, £(0), and 
a(t) = Bu(t’) is the forcing vector.’ Note that if A is a 
constant matrix, 


fe), (3) 


and (2) can be written 
E(t) = B(t)Eo +f @(t — t')a(t’)dt’. (4) 
0 


Because it is required that #(t) =&(f) at some future 
time, t=¢,, and that the time interval (0, ¢,) is to be a 
minimum, a requirement is that 


-1(1,)é(t,) — bo = f  e10)a)d’, (8) 


for some ®. Thus, the aim is to find the smallest ¢t, >0 for 
which an allowable control vector ® is, under the 
linear mapping given by (5), taken into the vector 
®-'(t,) E(tr) — Eo. 

The following two theorems concerning minimum- 
response-time control can be easily established for the 
process governed by the differential equation (1): 

Theorem I—lIf there exists a time optimal controller 
for the process (1), with #€Q, then there exists a time 
optimal relay controller, z.e., the control function can 
be on the boundary of Q for all ¢ in the time interval 
(0, ¢,). (A simple extension of the proof presented by 
Bellman, Glicksberg and Gross [23] will establish this 
theorem.) 

Theorem II—If &(0) is in the region for which there 
exists a time optimal controller, and if the ” character- 
istic roots of A, as seen in (1), are real and distinct, the 
components of the time optimal control vector @ 
should each change sign, at most, »—1 times in the 
time interval (0, ¢,). (A proof of this theorem is also out- 
lined by Bellman, Glicksberg and Gross [23 ].) 

A sufficient condition, if A is a constant matrix, for 
insuring the existence of a time optimal regulator is 


’ The matrix &() is the solution of the matrix system 
(0) = 7, 


where J is the unit matrix. If A is a constant matrix, the matrix ® 
is the matrix exponential 


@ = ct, 
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that the matrix A of (1) have characteristic roots with 
negative real parts. This is not a necessary condition, as 
is well known [25], for even in the case of unstable 
characteristic roots there may exist a region x over 
which it is possible to regulate optimally. An example of 
a situation where it is not possible to optimally regulate 
would be when S-1B=0. [See (10).] In this case no 
matter how much power is available or how intelligently 
used, the region x will be empty if any of the character- 
istic roots of the matrix A have a positive real part. 

Theorem I is very important in the synthesis of the 
minimum-time-response controller, since it enables the 
indicated integration of (5) to be carried out. The 
results of the integration are then used in the procedure 
for finding the optimal control function. The contribu- 
tion of this paper is to show how it is possible to go from 
Theorem I to a set of necessary equations from which 
the optimum control function can be found by means 
of a control computer. It will also be indicated how the 
switching equations, in terms of the state variables, can 
be found for the lower-order processes. 

Theorem II is not important to the method presented 
here, but if the conditions under which it holds are met, 
it will represent a major simplification of the method. 

In the work that follows it will be assumed that the 
matrices A and B are constant and that A has stable 
characteristic roots. Also it is assumed that £&(t)=0 
Neither of these assumptions is necessary in the use of 
the method presented. (The modifications for the more 
general case are obvious.) 


SINGLE- DEGREE-OF-FREEDOM PROCESSES 


The result presented in this section is actually a 
special case of the more general result presented later in 
this paper. From an engineering point of view, it appears 
that the simpler case should be studied first. 

The equation for the single-degree-of-freedom process 
can be written as 
d"x " hagiiage 


10 earn eo 
7 ae 


Wa + b,« = u, b; = constants, (6) 
fe 


which can be put in the form of (1) by letting 


= X 
and 
41 = Xe 
Xo = X3 
fle ole pe ee bh rieetae) (7) 
Thus 
dé 
5 = 48+ Bi 
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where 
0) i Oh 0 x4 
epee 0 Ty 1) de ae en 
ye eee ah, X 
0 0 0 0 
Re 0 0 0 
B= : : “=|. 
; 0 0 
0 0 1 U 


If distinct characteristic roots are assumed, the solu- 
tion to (6) can be written in terms of principal coordi- 
nates, v;, (v; is a linear combination of the original co- 
ordinates x;) as 


t 
a(t) = ety + J eM SHA(H) dr, (8) 


0 


where SS is Xm nonsingular matrix such that 
Gyn er Op raced) 
@ Gy Woo oO 
C= SAAS Sih tele 
; 0 a 


and the a; are the m distinct characteristic roots of the 
process equation. For optimum response it is required 
that 3(t) =0 when ¢t=t,, the response time; thus, 


ty 
Oia f erat at’, (9) 
0 


where for a single-degree-of-freedom process as previ- 
ously hypothesized, 


| Fe C12 eG “|: 
[ee (py 2 ONG) 


g(t) = Sot Serie 
R(t) ; als | 
E les Cro see G fp 
or 
2n Cnn 
nt) =e 
gill) = Ts] U, go(t) = Wists » &n(t) Rane (10) 


By Theorem I, | wi| =m; for all t in the interval (0, ¢,), 
and it will be assumed on the basis of Theorem II, that 
a finite but unknown number of control function switch- 
ings will include the time optimal controller. (In some 
cases it is possible to restrict the range of initial condi- 
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tions to require at most n—1 switchings.) The integra- 
tion of (9) can now be performed, as follows: 


te Pa 
Vm (0) = f emt’ g(t) dt’ = f emer" dt’ 
0 0 


to a 
— f Bre eam dt + ean + { &mer t dt! 
t re 


1 


= ott [1 — Deum + 2 eam te + tehet vs AE ZDeamtk ae eumtr|. 
am 
OES ie SETI EHS SO 
(Gs ey) et; ka > O° (11) 


The first question to be answered is: what is the smallest 
t,=0 for which the ” equations of (11) are satisfied? By 
answering this question, the present relay sign is auto- 
matically determined, 1.e., it will be the same as the sign 
preceding the first term, 2e*%, which has a nonzero f;. 
To get at the minimization of ¢,, one equation of (11) is 
solved for ¢,, in terms of the other switching times, ¢;, 


be) (12) 


and ¢, is eliminated from the remaining n—1 equations 
to obtain »—1 constraint equations, 


Go GAli, to,n4 oot) =. 0, pe lydia, 0 — 1. 


lie = F(t, to, : 


In addition we must have the inequalities 


ON ee ot pe ie (13) 


The minimizing procedure is then to look for the 
minimum of F(ty,- ++, t) over the set, as given by 
Gi(t1, tz, -- +, ¢), where the #,;’s are restricted by 
Opie a1) re ip. 920 Go this and not “be com: 
mitted as to the independent variable, a method devel- 
oped by Lagrange will be used. 

First form the function 


f= F(t, lo, at me) ae \iGi(t1, to, s 
+ Nannies to, ae ei’ tk), 
where the A; are constants, as yet undetermined in value. 


Treat h, f2, - - - , f, as independent variables and write 
down the necessary condition for an extremal, 


-- jth) te: 
(14) 


(15) 


subject to the constraints G;(é1, f2,--++, &) =0, 7=1, 
2,°°:+,m-—1, and where the ?;’s are restricted as fol- 
lows: 


OS iE GS OORT RS (16) 


It is now obvious how to obtain the proper relay sign. 
Simply mechanize (13) and (15) to compute the mini- 
mizing set of ¢;’s, and then select the present relay sign 
as previously indicated. These equations can be solved 
on-line, with either analog or digital computer type 
equipment [1]. Note that if Theorem II is applicable, 
the set of switching times to be considered has, at most, 
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n elements, 7.e., n=k. Note also that information in 
addition to the present relay sign will be obtained by 
this procedure, e.g., a prediction of when the next relay 
switch will occur. Finally, note that the procedure taken 
by the Russian authors [28]-[33], [35]-[38] and by 
Desoer [34] requires the information as presented here 
to obtain the optimal forcing function; 7.e., if given the 
initial conditions of the adjoint equation, they can ob- 
tain a function which has zeros at the switching times, 
or vice versa. Their equations can be solved for either 
the switching times or the initial conditions of the ad- 
joint equation, but not for both. 


GENERAL LINEAR PROCESS 


The general linear process is governed by (1), which 
has a solution as given by (2). Consider the necessary 
equation (5) which must occur at t=?,, the response 
time. 

Call the elements of the matrix ®~'(t), aij, 7.e., 


ais(t) aye(t) + + + ain(t) 
€ ee oa (!) atoe(t) > ~+ aen(t) 
(17) 
ani(t) Ano(t) eat Ann(t) 
and write 
bist + Dine = 4 Oi thn 
boty + boot, + + + + + Donttm 
a(t) - 721 1 20 2 (18) 
baits SF byte a © SPs SF Dnmtbm. 


where each | u;| =m; for all t in the interval (0, ¢,). The 
integration of (5) is now performed in parts, assuming 
k; switching times for each independent forcing func- 
tion, 2;. 


—2(0) = ii Bull yua() dt! +5 if “Bia(l")ao(t") de! + or 
0 0 
+f ‘Bim(t’) u(t!) dl’, 
0 
—2,(0) = i ‘Bolt! an(t!)dt + f “Boat! us(")dl rn 
0 0 


a f (oe (¢’) Um(t' dt’, 
0 


OE (he f ‘Bart’ (t)d0! + i Bat P ear eae 


ty 
+f Burt!) Ul Oat. (19) 
0 
Carry out the integration just for the first equation, let- 
ting ¢11, tia, fis, - > + , tm, be the ky switching times for 


Uys botydaas © <5 bony, © 2 ey ee 
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=v (0). = = Bu (t’)dt’ sf Bill’ dt’ + . Jee] 
0 


oils 


tm. 
si i Bin(t’) dt! ty 
0 


with similar equations for x2(0) - - - x,(0). Now solve 
one of the above equations, as before, for ¢,: 


te = F(tu, tie >> + ti, ter - Mat, (240) 


and eliminate ¢, from the remaining »—1 equations to 


obtain —1 constraint equations, 
G; = Gi(tu, nyo 2 oak) — th I. OG dig ae if (22) 


To find the minimizing set of switching times, pro- 
ceed as before, forming the function. 


f= FHAGi + d2Go+ +--+ An1Gru, (23) 
from which we calculate 
of et VAR 
a= =O; = 1,2,---, km’, (24) 
Ot;; ia V2 ><, Mm 
subject to the constraints G;, 1=1, 2:---mn—1 and 


where the f;;’s are restricted as follows: 


Osta Ss fs ee 1= 1, 2, “5M. 


The sign of each relay is then chosen to correspond to 
the appropriate term involving the first nonvanishing 
t;;, For more general results and the use of other opti- 
mum criteria see [41]. 


A SECOND-ORDER EXAMPLE 
Consider a particular second-order process as de- 
scribed by the differential equation, 
41 = AyoX%2 + bum 


(25) 


Lo = dex + A22%X x2 boots. 


Assume that this equation has real distinct characteris- 
tic roots —a, —b. The fundamental matrix ®(¢) is then 


J be = ae ° eat fad et | 
b—a b— a 
—abe-** + abe** ~—ae + be "| 

b-—a b—a 


Because A was assumed to be a constant matrix. 


[ bet — ae gtr et ] 

b—a ba 
anil — _ => iY (27) 
Se aie + aber! —aert + | 


Ca b—a 


Bio(t')dt’ — “ntl )dt’ + - 
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lift 
te 
Brat) | 

lore 

tme tr 
Bim(t’) dt! + pial '4 i mata |, (20) 

t 


t 
mi bm k 
mn mk m! 


Writing out (5) for the regulator problem, 


tr ber’ ek ae’ eat aE ebt! 
¥1(0) => f{ | am by, + ‘ia ah be | are 
0 b—a b—a 
TT —abe*” + aber!’ 
%2(0) = i | bitty 
0 b — a 
—aer"’ + beet’ 


bu | di. (28) 
b=@ 


which can, in turn, be written in principal coordinates 


tr 


€°! (abyy M4 + boots) dt’, 


e 
= 


= ax,(0) + x(0) ={ 


0 


ty 


7! (bbi 14 + boots) dt’. 


l| 


2 bx,(0) + x2(0) -{ 


0 


(29) 


Because each | 20 =1 (letting the matrix B carry the 
magnitude) for all ¢ in the interval (0, t,), and the char- 
acteristic roots —a, —b are assumed real and distinct, 
the above integration can be carried out in n—1 parts, 
assuming an initial sign for the forcing terms, #, and 2». 


—aby, 


b 
[1 — 2eb 4 ebtr] — = [1 — er + ert], 


—bb b 
11 [1 — Qe 4 ert] asa v2 [1 — Jette + evtr| 
a a 


2bb4, 2b» 
—_.—| —— evi — oe eta 
a a 


(30) 


vg = 


Solving for ¢, 


[=| |= 
—v,— | —— ebtt pee ae ete 
| ihe : : (31) 
t;=— log) — : 
eee 
b b 


Call the first equation of (31), F(t1, t2), and eliminate 1, 
from the second equation to find the constraint equa- 
tion 
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G(ti, t2) 
\—"| || 1/a 
—V_. — see ewth ~) pee es e%t2 
0) 1+ ; 2 
- bb, bas 
—+ 
a a 
| [= 1/b 
—Vy — ebtt ebte 
b 
ms yeas cole 
aby, bo» 
+ 
b b 


alg ios eae (32) 


The minimizing procedure is then to look for the mini- 
mum of F(t, t2) over the set in the first quadrant as 
given by G(éi, t:). Let 


f = F(hy, te) + AG(A, te) 


1 
= — log [7] + A[zZ]}*/¢ — a[zz]}2”, (33) 
a 
Thus, 
— aty 
of ss 2bbi1e | 1 2 ate | 
84 = bb + B22 L[Z] 
—2abi1 
_ ize] Jor = 0, (34a) 
aby Be boo 


af [ 1 
eta ees 2b to x Jé (1/a)—1 
als 22€ [7] = [7] 
—2boo 
— Vj [ZT] /e)-1 |= | Star=an()) 34b 
(Leroy [Jom = 0, Gan) 


Eqs. (32) and (34) can be used to solve for the switching 
times and, in effect, the proper relay signs by means of 
a control computer. 

To elucidate further the importance of the procedure 
presented, consider now the problem of determining the 
switching boundaries for the two forcing functions in 
the particular second-order example being considered. A 
simple algebraic simplification of (34) leads to the equa- 
tion expressing the relationship between the switching 
times which occur at the minimum. 


aq ]1/(a-2) 
ty; — te = lo =| 3 
1 2 = b 


Thus, in all cases, the switching times ¢; and f2 differ by 
at most a constant. To plot a curve for the switchings of 
u; and us, select a=1, b=2. First note, in this case, that 
ty is larger than ft. by, at most, log 2, and that if the aim 
is to find the switching boundary along which both 
forcing terms are positive, there will be no switches in 
the forcing terms, 17.¢., 4: =t2=0. The second condition 


(35) 


(36) 
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A simple test of the direction of movement of a point on 
this boundary will determine which portion is the re- 
quired switching boundary; 1.e., it should be the portion 
of the curve along which the point moves toward the 
null point. 

By similar reasoning, the other three boundaries can 
be determined. Condition 1 has not yet been imposed 
for the minimum which, for this example, is given by 
(35) with a=1, b=2. It has been noted before that 4, 
is larger than fy; thus, to find the switching boundary 
for us, it is just necessary to consider the locus of points 
that are log 2 units of time away from a switching curve 
as determined previously. To find the switching equa- 
tion for uw, consider (29), which can now provide a rela- 
tionship between the optimum switching times f; and fg 
as given by (35). The result of this substitution leads to 
the equations, 


Y= 


1 bos 
bar 


v2 = — 2bn[—3 + e*] — boo[—1 4+ eX]. (37) 
Eliminating ¢,, the response time, obtain the following 


equation which is part of the wz: switching boundary: 


Ee + doe 201 | 
bir + doe bir + doe 


6b11 + doe V2 | 
= + - (38) 

2611 Se boo 2611 me bo» 
Fig. 1 indicates, for the particular second-order example 
considered, what the complete switching boundaries for 


uz and u2 look like in the phase plane composed of the 
principal coordinates v; and 2. 


y, | 20 


Ka “Ky 


-20 


Fig. 1—Optimum switching curves for second order process with two 
independent forcing terms, bi =2, b22=5, a=1, b=2. 
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CONCLUSIONS 


A method has been presented which leads in a system- 
atic manner to the desired time optimal forcing func- 
tion. The method lacks esthetic beauty and does not 
give any general insight as to the over-all structural 
properties of the problem, but it does lead to the re- 
quired engineering solution, which is to know the pres- 
ent sign of each relay forcing term. Admittedly, the 
method requires a great deal of care in computing the 
solution to the transcendental equations obtained under 
the constraints imposed; this difficulty is correctly 
amplified in the discussion which follows this paper. 
However, as [1] shows, it is possible with presently 
available computers to use the method presented here 
to find, in real time, the time optimal control function 
for a large class of present-day high-order dynamic proc- 
esses. 
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Discussion 


A. M. Hopkin (University of California, Berkeley): 
The ideas expressed in this paper are stated in more ele- 
gant and general terms, but bear some resemblance to 
previous work done by the discussor.! On the basis of 
this experience, the following comments seem pertinent. 

In (13) of the paper two kinds of constraints are speci- 
fied, one set up in the form G;=0, and the other in the 
form O<fi Stem = lp St Phe: first constraint is 
taken care of by (14) (Lagrange’s multipliers), but the 


1 A. M. Hopkin and P. K. C. Wang, “A relay-type feedback con. 
trol system designed for random inputs,” Trans. AIEE, vol. 59. 
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second is not. If a computer is programmed to minimize 
the Lagrangian function of (14) by setting derivatives 
with respect to various f; variables equal to zero, it is 
possible that a solution might be obtained which did 
not conform to the second constraint. It is also con- 
ceivably possible that for a particular f¢; and set of initial 
conditions there would be no zero of the partial deriva- 
tive of the Lagrangian function with respect to ¢; 
within the constrained range of ¢;, yet the function could 
be minimum at one end of that range. It is felt that these 
factors require careful consideration before any control 
computer is built. 

Another possibility which exists is that the Lagrang- 
ian function as seen in the multidimensional ft; space 
might be a surface with many humps and hollows, so 
that the minimization equations (15) might have many 
solutions. This event would be likely if some of the char- 
acteristic roots of the system equation occurred in com- 
plex conjugate pairs. 

In this case, the computer would have to search 
through all possible solution sets to find the particular 


IRE TRANSACTIONS ON AUTOMATIC CONTROL 


September 


one which met the constraints and gave a true minimum 
ty. Since the computer must solve simultaneous tran- 
scendental equations, the total solution might be quite 
slow in terms of the time scale of the controlled process 
for any economically designed computer. 

The concept of designing a control system on the 
basis of minimization of the time required to bring the 
controlled process from an initial state to some desired 
state was first introduced almost ten years ago. Contin- 
uing contributions are being made to the general theory. 
Unfortunately, it seems that the details involved in the 
design of actual systems utilizing this theory are 
troublesome, so that very few practical systems have as 
yet been built. The main use of the theory so far has 
been to provide a model for comparison. 

Author’s Comments: I would like to thank Professor 
Hopkin for his discussion of my paper. His discussion 
correctly amplifies a difficulty which is encountered in 
solving the time optimal control problem. In part, these 
difficulties have been overcome; the solutions so far ob- 
tained are reported by F. Smith [1]. 


Statistical Design of Digital Control Systems” 
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Summary—This paper is concerned with the problem of sta- 
tistical design of digital control systems. A simplified design pro- 
cedure is introduced. The minimum mean-square error is adopted 
as the performance criterion; and the system design is based upon 
the Wiener-Kolmogoroff theory of optimum filtering and prediction. 
The treatment of this optimization problem makes use of the modi- 
fied z-transform technique, and is presented in close parallel with 
the procedure for the statistical design of linear continuous-data 
control systems. By making use of some important properties of the 
auto-correlation sequence, the cross-correlation sequence, the pulse 
spectral density and the pulse cross-spectral density, the statistical 
design of digital control systems can be carried out in a simple and 
systematic manner. 
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INTRODUCTION 
D URING the past decade, the analysis and design 


of digital control systems have been extensively 
studied. In the past few years, several papers'> 
have been written on the subject of statistical design of 


'G. Franklin, “Linear filtering of sampled data,” 1955 IRE Na- 
TIONAL CONVENTION RECORD, pt. 4, pp. 119-128. 

2 R. M, Stewart, “Statistical design and evaluation of filters for 
the restoration of sampled data,” Proc. IRE, vol. 44, pp. 253-257; 
February, 1956. 

* A. S. Robinson, “The synthesis of computer-limited sampled- 
data control systems,” ATEE Special Publication T-101, pp. 77-78; 
May, 1958. 

*S. S. L. Chang, “Statistical design theory for digital-controlled 
continuous systems,” Trans. AIEE, vol. 77, pt. 2, pp. 191-201; 
September, 1958. 

5 J. T. Tou, “Digital and Sampled-Data Control Systems,” Mc- 
Graw-Hill Book Co., Inc., New York, N. Y., Ch. 10; 1959, 
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sampled-data systems. The procedures proposed in 
these papers appear to be either involved or incomplete. 
This paper attempts to study the statistical design 
problem by an alternate approach so that a simpler de- 
sign procedure may result. Compact design formulas 
are derived in close parallel with the statistical design 
of continuous-data systems.*-§ Optimum pulse-trans- 
fer functions are determined on the basis of the mini- 
mization of both the mean-square error and the mean- 
square sampled error. The design procedures are dis- 
cussed along with physical realizability and system con- 
straints. In this paper the input signals and the noise of 
the system are assumed to be stationary random func- 
tions. 


SOME PROPERTIES OF CORRELATION SEQUENCE 
AND PULSE SPECTRAL DENSITY 


When the minimization of the mean-square error is 
adopted as the performance criterion, the design of dis- 
crete-data systems is usually initiated with the descrip- 
tion of the signals by the correlation sequences and 
pulse-spectral densities.® The definitions of auto-correla- 
tion sequence, cross-correlation sequence, pulse-spec- 
tral density, and cross pulse-spectral density are given in 
the literature. It has been shown’ that the correlation 
sequence and pulse-spectral density possess the following 
properties: 

1) If the response of the pulsed-data system G(z) to 
an input x*(t) is y*(t), then the response of this system 
to an input ¢,:(R7) is ¢:,(kT), and the response of this 
system to an input-¢,,(R7) is d,,(RkT), where $,2(kT) 
and ¢,,(k7) denote the autocorrelation sequence of the 
pulsed signals x*(t) and y*(t), respectively, and ¢z,(kT) 
and ¢,:(k7) are the cross-correlation sequences between 
these two pulsed signals. 

2) The pulse-spectral densities for the input and the 
output sequences of a pulsed-data system G(z) are re- 
lated by 


dyy(2) ae G(z)G(z“!)dz2(2), (1) 


where ¢,,(z) is the pulse-spectral density for the input 
sequence x*(t), and ¢,,(z) is that for the output sequence 
y*(t). 

3) Assume that two pulsed-data systems Gi(z) and 
G.(z) are subjected to input sequences x;*(t) and x»*(t), 
respectively, and that the output sequences in response 
to these inputs are y:*(t) and y2*(t), respectively. Then 
the output sequence of the system Gi(z) is $x,,(kT) 
when it is subjected to an input sequence $2,2,(k7), and 


6H. S. Tsien, “Engineering Cybernetics,” McGraw-Hill Book 
Co., Inc., New York, N. Y.; 1954. ; 

7]. H. Laning and R. H. Battin, “Random Processes in Auto- 
matic Control,” McGraw-Hill Book Co., Inc., New York, N. Y.; 

6. 
ae G. C. Newton, L. A. Gould, and J. F. Kaiser, “Analytical Design 
of Linear Feedback Control,” John Wiley and Sons, Inc., New York, 
ie Mee TOkye ae 
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the output sequence of the system G2(z) is $y,,,(RT) 

when it is subjected to an input sequence $y,7,(R7). 
4) The cross pulse-spectral densities for the input 

and the output sequences of two pulsed-data systems 


Gi(z) and G:(z) are related by 
by ua(2) = Gi(271)Go(2) bey29(Z), (2) 


where ¢2,2,(g) is the cross pulse-spectral density for the 
input sequences x,*(¢) and x2*(t), and ¢y,,.(z) is that for 
the output sequences y,*(t) and y»2*(Z). 

Furthermore, it can be shown that the following in- 
put-output relationships also hold: 

5) If the response of the linear system G(s) to an 
input sequence x*(¢) is y(t), then the “modified” pulse- 
spectral density for the output y(¢) is related to the 
pulse-spectral density for the input x*(t) by 


by (2, m) = G(z, m)G(z-", m)dbzz(z), (3) 


where G(z, m) is the modified z transform associated 
with G(s). 

6) Assume that two linear systems G(s) and Go(s) 
are subjected to input sequences x1*(t) and xe*(#), re- 
spectively, and that the system outputs in response to 
these inputs are yi(¢) and y2(t) respectively. Then the 
“modified” cross pulse-spectral density for the output 
signals yi(¢) and y(t) is given by 


byiya(2; m) rz Gis, m)Go(z, Mm) bz25(2) (4) 


where G,(z-!, m) and G2(z, m) are the modified z trans- 
forms associated with Gi(—s) and G.(s), respectively, 
and @$2,7.(2) is cross pulse-spectral density for input se- 
quences x;*(f) and xo*(é). 

7) If the response of linear system G(s) to an input 
sequence x;*(t) is y:(t), and the response of linear sys- 
tem G.,(s) to an input signal x2(t) is ye(t), then the 
“modified” cross pulse-spectral density for the output 
signals yi(t) and yo(t) is given by 


bywolZ; m) ee GiGi(z, M) bx29(2) 5 (5) 


where G:G.(z, m) is the modified z transform associated 
with Gi(—s)G2(s). 

The above properties of correlation sequences and 
pulse-spectral densities will find much use in the opti- 
mum design of digital control systems based upon the 
Wiener-Kolmogoroff theory. 


Optimum PuLSED-DATA COMPENSATOR FOR MINIMUM 
MEAN-SQUARE SAMPLED ERROR 


Thestatistical design problemcan be readily visualized 
by referring to Fig. 1. Go(z) is the over-all pulse-transfer 
function of a digital feedback control system which is 
to be optimized on the basis of minimum mean-square 
sampled error. The sampled input and output of the 
system are assumed to be 


r(nT) = 7.(nT) + 1r,(nT) (6) 
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Fig, 1—Digital control system with sampled output. 


and 


c(nT) = ¢.(nT) + cr(nT), (7) 
respectively, where the subscript s refers to the signal 
and the subscript m refers to the random noise. Also 
assume that the desired pulse-transfer function be 
Ga(z) and the desired output sequence be ca(nJ7). The 
error sequence of the system is then given by 
e(nT) = c(nT) — ca(nT) 
= ¢,(nT) + Cr(nT) — ca(nT) (8) 


and the mean-square sampled error is 


e?(nT) = cs?(nT) + cn?(nT) + ca2(nT) + co(nT) en(nT) 
+ ¢n(nT)c.(nT) — Cn(mT)ca(nT) — ca(nT)cn(nT) 


— ¢.(nT)ca(nT) = ca(nT)c.(nT). (9) 
In view of the relationships 
ah Zs" 1 
x?(nT) = $22(0) = — ¢ bza(2)2~*dz (10) 
21] T 


and 


Poi Pee ed ED il 
AT IVAT) = bn(0) = 5 ¢ eu (Oya pe edi 
TT iM 


where the contour of integration is the unit circle, the 
mean-square sampled error given in (9) may be written 
as 


See) | 
e?(nT) = ae [bege,(2) = Were, (2) = Peaca(Z) am Pegcn (2) 
Fr 


+ Penes(Z) — PencaZ) — Degen (2) — be,ca(2) 

— Hege,(2) |z~dz. (12) 
By making use of the properties of the correlation se- 
quences and the pulse-spectral densities given in (1) and 


(2), the above equation for the mean-square sampled 
error may be reduced to 


(nT) = wif ours (13) 


in which 
ee(2) = [Go(z) — Ga(z)|[Go(z) — Ga(z) ]o,,-,(2) 
+ Go(z)[Go(z-!) — Gaz?) ]bryrn (2) 
+ Go(2-) [Go(z) — Gal) ]brar,(2) 


+ Go(Z)Go(2-") brnrn (2) (14) 


AUTOMATIC CONTROL Septem ber 


With reference to Fig. 1, it is seen that the over-all 
pulse-transfer function is related to the pulse-transfer 
function of the controlled system G(s) by 

Go(z) = W(z)G(z) (15) 
where 
D(z) 
d 
1 — D(z)G(z) 


W(z) = (16) 


and D(z) is the pulse-transfer function of the pulsed- 
data compensator. Substituting (15) into (14) yields 


Peel) = [W(2)G(2) — Galz) [We GE) — Gale) ]$r41.(2) 
+ W(z)G(z) [We ) Ge) — Galz) ]Orarn(2) 
+ We )G(e) [W(@)G@) — Gal2) |orar(2) 
+ W(e)W (2G) G2) bran (2) + (17) 


The design objective is the determination of the pulse- 
transfer functions W(z) and D(z) in such a way that the 
mean-square sampled error given by (13) is reduced to 
a minimum. It follows from (16) that D(z) is given by 


W(z) 
1 — W(2)G(z) 


D(z) (18) 


The minimization problem can readily be solved by 
applying the calculus of variations to the integral of 
(13) with ¢..(g) given in (17). To determine the condi- 
tion for minimizing the mean-square sampled error, the 
pulse-transfer function W(z) is given a small variation 
An(z), and W(zg7') is given a small variation \n(z7'). As 
a result the mean-square sampled error e?(n7) under- 
goes a variation 6e2(n7). n(z) is an arbitrary function of 
z, and )d is a parameter. Replacing W(z) by its neighbor- 
ing function W(z) +An(z), and W(z) by its neighboring 
function W/(z~!)+An(z7!), the mean-square sampled 
error becomes e?(n7) +6e2(n7T). Then e2(n7) is a mini- 
mum, if 


d[e?(nT) + 6e(nT)| 
dn 


ay. (19) 


A\=0 


Following the operations involved in (19) and simpli- 
fying leads to 


= $ eG) (W(2)G(2)$or(2) 
— Gal2) [brar4(2) + rare) ]} o-¥de 
+ af 260 { W(z)G(z) bn(2) 


— Gaz) [brr(2) + brinn(2) |}2dz = 0 (20) 


where 


brr(Z) = Prarg(Z) + Prarn(Z) + Prern(Z) + brar,(Z). (21) 
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To determine the optimum W(s), the z transform ¢;;(z) 
is factored into two parts, one part having zeros and 
poles inside the unit circle and the other part having 
zeros and poles outside the unit circle. Thus : 


ber() = burt(2)bn-(2). (22) 
The zg transform G(z)G(z~) is similarly factored: 
G2)GE) = [G@GE)FIE@GE)/ — (23) 


where the first factor contains zeros and poles inside the 
unit circle and the second factor contains zeros and 
poles outside the unit circle. Then (20) may be written 
as 


1 
en (2) ber) [E@Ee 1 [G2)GE) FW 2) bat (2) 


21] 


as G(z)Ga(z) [by,.(2) + =) 12h se 
[G(z)G(z) |- be (2) ; 


1 
“3 ami g ale) ber (@) [E@Ge)|* {Ie@oeyee.-@ 
_ GO)GAEN Fr) + orm ON 17 _ 9 os 
[6G@ce)}e+@ J 


For physical realizability, W(z) and n(z) can have no 
pole outside the unit circle of the z plane, and W(z"') 
and 7(z~) can have no pole inside the unit circle. Thus, 
the function 


1) bre (2) [G@)G(e) > 
contains poles outside the unit circle only, and the func- 
tion 
n(2) br (2) [G(z)G(z-}) |* 


contains poles inside the unit circle only. Consequently, 
the above equation may be further reduced to 


a nz) bre (2) [G(2)G(z-) |r | [G(z)G(z-) |+W (2) br (2) 
7] T 


is (see + ae | 
[G(z)G(z) |b, (2) + 


+= n(dorA COG) [6G We) 
Ti] Tr 
_ {eoen [Sraro(2) + waa be pees 
[G@)GE) |*b-r*@) s 


in which the symbol { |, designates the operation of 
picking the part of a function of zg with poles inside the 
unit circle, and the symbol { |_ designates the opera- 
tion of picking the part with poles outside the unit circle. 
For instance, in F(z) = {F (z) tet { F(z) pat the terms 
with poles inside the unit circle are lumped in the (+) 
part, and the terms with poles outside the unit circle 
are collected in the (—) part. 

~ Now, if W(z) is the optimum pulse-transfer function 


ys 


Tou: Statistical Design of Digital Control Systems Z9e 


for minimum mean-square sampled error, (25) must be 
satisfied for arbitrary n(z). Thus the optimum W(z) is 
given by 

(G(2)Galz) [brors(2) + Pras) J) 


L [e@e@oFe-@ 54 _ 
[G()G) |+bn*(2) 


W(z) =- 


(26) 


Once the pulse-transfer function W(z) is determined, the 
required pulse-transfer function of the compensator 
D(z) follows immediately from (18). If the pulse-trans- 
fer function G(z) of the controlled system contains no 
zero and no pole outside the unit circle of the z plane, 


[G(2)G(z) |* = GG) (27) 
and 


[G@GE))> = Ge"). (28) 


Consequently, under such conditions, (26) can be 
simplified to 


Ga(z) [brsrs(2) a ell 
Fe 


ie drr (2) 
W(z) = Gwe. (29) 


Eqs. (26) and (29) form the working formulas, from 
which the optimum W(z) can be readily determined. 
The determination of W(z) and D(z) simply involves 
the evaluation of zg transforms and simple algebraic 
manipulations. Clearly, with the aid of these two equa- 
tions, the design of digital control systems for minimum 
mean-square sampled error can be carried out syste- 
matically in a simple manner. 


Optimum PuLSED-DATA COMPENSATOR FOR 
Mrnimum MEAN-SQUARE ERROR 


Now, a slightly different problem is to be considered. 
This section will present a procedure for the design of 
the pulsed-data compensator to minimize the mean- 
square error of the system. To determine the optimum 
pulse-transfer function in a way as simple and system- 
atic as the design procedure presented above, use is 
made of the modified z-transform technique. 

With reference to Fig. 2, Go(z, m) is the over-all 
modified pulse-transfer function of the digital feed- 
back control system which is to be optimized on the 
basis of minimum mean-square error. The input and 
the output of the control system are assumed to be 
r(t) =rs(t) +rn(t) and c(t) =c,(t) +en(¢), respectively. Lt 
is also assumed that the desired system transfer func- 
tion is Ga(s) and the desired output signal is ca(t). The 
system error is then given by 


e(t) 


I 


c(t) — calt) 
c(t) + en(t) — ca(t). (30) 
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Fig. 2—Digital control system with continuous output. 


In terms of samples, 
e(nT, m) = c.(nT, m) + cr(nT, m) — ca(nT,m). (31) 


It follows from (31) that the mean-square sampled-error 
for a particular value of m is given by 


e?(nT, m) 


= ¢.(nT, m) + Gr2(nT, m) + ca?(nT, m) 
+ (nT, m)cn(nT, m) + cr(nT, m)c.(nT, m) 
— ¢.(nT, m)ca(nT, m) — ca(nT, m)c.(nT, m) 


— €n(nT, m)ca(nT, m) — ca(nT, m)cn(nT,m). (32) 


Following the argument of the previous section, (32) 
may be reduced to 


PUT) = 5 fdas, mers, (33) 
2ajJ/ r 
in which 
Pee(Z, Mm) 
= degeg(Z,M) + Penen(Z, mM) + Pegca(Z, Mm) + He,e,(z, m) 
F Pencg(2, ™) — Pegea(Z, ™) — degc,(Z, m) 
— Penea(2, ™) — Hegen(Z, M2). (34) 


In view of the properties of the pulse-spectral densities 
given in (3) and (5), @ce(z, m) may be expressed in terms 
of the modified z transforms of the system transfer func- 
tions. Thus, 


bee(z, m) = [Go(z, m)Go(z-1, m) — GoGalz, m) — Gian m) 


G,(z, m) 
ao 
+ [Go(z, m)Go(z-, m) — GoGulz, m) |brer,, (2) 
+ [Go(z, m)Go(z, m) — GoGalz, m) |]bryr,(2) 
+ Go(z, m)Go(2-!, M) brary (2) (35) 


where GpGa(z, m), GoGa(z, m) and G,(z, m) are the 

modified z transforms associated with Go(=s)Ga(s), 

Go(s)Ga(—s) and Ga(s)Ga(—s),,7,(s), respectively. 
Examination of Fig. 2 reveals that 


Go(s) = W*(s)G(s) 


(36) 
where 

D*(s) 
1+ D*(s)G*(s)_ 


W*(sj= (37) 
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Taking the modified z transform of both sides of (36) 
yields 


Go(z, m) = W(z)G(z, m). (38) 


By substitution of (36) and (38), ¢-<(z, m) given in (35) 

becomes 

Pee(2, Mm) 

= [W(2)W(z)G(z, m)G(z-1, m) — W(z)GGi(z, m) 
— W(z)GGi(z, m) + G,(z, m) /brer,(2) |brerg(Z) 
+ [W(e)W(e)G(z,m)G(=-1,m) — W(2)GGa(z, m) |brora(2) 
+ [W(2)W(z)G(z, m)G(z, m) — W(z-)GGalz, m) |brar,(2) 
+ W(z)W(z)G(z, m)G(z-!, m) brar, (2). (39) 
It has been shown that the mean-square error e?(#) can 


be determined from the mean-square error sequence 
e(nT, m) through integration;> that is, 


e-(t) = [ e0r, m)dm. (40) 


Combining (40) and (33) gives the mean-square error as 


1 


er vids f Peels, m)dm (41) 


on 


where 


1 
f hee(Z, mdm 
0 


= [W(2)W(e)Kols) — W(2)Ko(2) — W(s) K(z") 
+ K3(2)/brer4(2) ]brara(2) 
W (2) K2(2) |brpr,, (2) 


W (z“1) K» (z~}) ]brnr, (2) 


+ [W(2)W(e) Ko(z) — 
+ [W(2)W(e) Ko(z) — 


+ WWE) Kol) dre) (42) 
in which 
Kole) = ff Gls, mG, mam (43) 
K2(z) = if CCuls, m)dm (44a) 
Kile) = I CGila, mam (44b) 
Ks(@) = f° Gals, mdm (45) 


Since (41) bears a close resemblance to (13) of the 
preceding section, the minimization problem may be 
solved in a similar fashion. 

Application of the calculus of variations to the inte- 
gral of (41) yields the optimum pulse-transfer function 


je [Prare(2Z) + dears (2)] | 
Ko (2) Por (z) EM 
Kot (2) brr*(z) 


W(z) = (46) 
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Clearly, the pulse-transfer function of the required 
compensator D(z) follows immediately from  substi- 
tuting (46) into (18). Eq. (46) can serve as the working 
formula for the statistical design of digital control sys- 
tems based upon the minimum mean-square-error cri- 
terion. Since (46) involves simply the evaluation of ¢ 
transforms and modified z transforms, ordinary integra- 
tion and simple algebraic manipulations, the determi- 
nation of the pulse-transfer function of the required 
compensator D(z) is indeed a simple matter. 


ILLUSTRATIVE EXAMPLES 


Example 1: Design of digital compensator for mini- 
mizing mean-square error sequence. 
Referring to Fig. 1, assume that 


hee 
s?(s + 1) 
Gas)-= 1 T = 0.2 second 
1 
Grr, (®) = Angee dra, (@) = 0.1 


dramn(®) = Prar(w) = 0. 


Determine the transfer function of the digital compen- 
sator required for minimizing the mean-square error 
sequence. 

The design is initiated with the determination of the 
z transform associated with G(s), ¢,,,,(5), @r,r,(s), and 
¢,r(s). It is found that 


0.019(z + 0.895) 


eRe Hie 30.819) 


— 0.20132 
(z — 0.819)(z — 1.221) 
Prar,(Z) = 0.1 
0.1(z — 0.261)(z — 3.791) 
(z — 0.819)(z — 1.221) 


Prar,(Z) = 


Prr (z) oF 


From the above equations, one obtains 


0.316(z — 0.261) 


OME 819) 
Oe Waste 3790) 
$e) ee i951) 
and 
ie Chee anode 
eae 7 pBO819 


Making use of (26), 


29.25(z — 1)(z — 0.819) 


We) = 9.261) (@ + 0.895) 
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Substituting into (18) yields the transfer function of the 
required compensator as 
29.25(2? — 1.8192 + 0.819) 


D(z) : 
(2? — 0.082 — 0.73) 


Example 2: Design of digital compensator for mini- 
mizing the mean-square error. 
Consider the system shown in Fig. 2. Assume that 


iheop see 
s(s + 1) 
Gas a1 i = 0.2 second 
bre (W) = Ona O01 
1 + w? 


Pr gry, (@) = Prprg() == (@), 


Determine the transfer function of the digital compen- 
sator required for minimizing the mean-square error. 

To design the system by following the above proce- 
dures, the z transforms Ko(z) and K2(z*) are first to be 
determined. Since the modified z transform associated 
with G(s) is 


1 (a= ies 2 
z 2(z — 0.819) 


G(z, m) = 


K((z) = Hf Cl, m)G(z—!, m)dm 


—0,0069(z — 0.305) (z + 3.274) 
‘ (z — 0.819)(z — 1.22) 


and 


K2(z—) 


I 


1 
if G(z-!, m)dm 
0 


—0.1072(z + 1.07) 
(g — 1.22) 


Factoring Ko(z) yields 


0.083(z + 0.305) 


OG ae ee aan atts 


—0.083(z + 3.274) 


ua (g — 1.22) 


Then 
0.0871 


et -_— 


Ko (2) rr (2) 
It follows from (46) that 


SO2(g — 0.819) 


Ws (= 0.261) (2 40.305) 
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Therefore, the transfer function of the required compen- 
sator is given by 
3.32(z — 0.819) 


D(z) = : . 
(z? + 0.0442 — 0.669) 


CONCLUSIONS 


This paper introduces a simplified procedure for sta- 
tistical design of digital control systems. The design 
techniques presented above are based upon two per- 
formance criteria, the minimization of mean-square 
sampled error and the minimization of mean-square 
error. These two criteria apply to purely digital control 
systems and digital-controlled continuous-data systems, 
respectively. The input signals and noise are assumed 
to be stationary random functions. In this paper the 
optimization problem is treated systematically by use 
of the modified z-transform technique. To simplify the 
design procedure, use is made of some important proper- 
ties of correlation sequence and pulse-spectral density, 
which are summarized above. Based upon the above 
two performance criteria, general design formulas are 
derived, which can readily be applied to the determina- 
tion of the pulse-transfer functions of optimum digital 
compensators. In deriving these design formulas, special 
attention is paid to physical realizability, nonminimum- 
phase type transfer functions and system constraints. 
Following the procedures presented above, the statis- 
tical design of digital control systems is made no more 
difficult than the statistical design of conventional con- 
tinuous-data control systems. Discussions reported in 
this paper are concerned only with digital control sys- 
tems having noise occurring in the input channel. The 
optimum design of digital control systems with noise or 
disturbance occurring at other points of the control loop 
will be discussed in a forthcoming paper. 


Discussion 


D. P. Lindorff, University of Connecticut, Storrs, 
Conn. The author has presented a clear and detailed 
derivation of the optimum pulsed filter in the complex- 
frequency domain, imposing various constraints. It 
may prove of some interest to point out that the results 
obtained can be derived directly by using the Bode- 
Shannon approach.! 

Thus in the model shown in Fig. 3, which is equiva- 
lent to that of Fig. 1, the optimum filter without regard 
to realizability is seen to be 


Ga (z) Orars (z) i 
rr (2) 


In this discussion ¢,.,;,,(z) =0 for simplicity. 


pGoe) = (47) 


1H. W. Bode and C. E. Shannon, “A simplified derivation of 
linear least-square smoothing and prediction theory,” Proc. IRE, 
vol. 38, pp. 423-424; April, 1950. 
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The optimum realizable filter is then determined 
directly to be 


1 ad 
Go + = ) 
on (t) | all bre (2) le 


where the term represented by the first factor is chosen 
to produce white noise. The term represented by the 
second factor is then taken to be the realizable com- 
ponent of the filter, which would produce the over-all 
optimum filter without regard to realizability. 


Gr by 
git Saal 
Ovnlod + 
e (nt) 


(48) 


If constraints are imposed upon G(z), resulting in 
(26), wherein Go(z) = W(z)G(z), the desired expression 
for W(z) can be obtained by the Bode-Shannon ap- 
proach, if the diagram of Fig. 3 is redrawn in the form 
shown in Fig. 4. Considering ¢,,,,/(z) and @,,,,/(z) as in- 
puts, the realizable optimum expression for W(z) is seen 


to be 
1 Gal2) One. (2) 
rats [ate 9 
dre (2) IL G(z) br’ ~(z) 
where it is understood that Ga(z)/G(z) replaces Ga(z) in 


Fig. 3. By using the following relationships, which are 
evident from Fig. 4, 


rere (2) = Pryr,(2)G(z)G(z-}) 
Prutn (2) = Prnrm(2)G(2)G(z-!) 
brr (2) = brr(z)G(z)G(z-), 


(49) 


(49) becomes 
Wopt (z) 


1 Ga(z)G(2) br 4r,(2) 
= —_ 50 
| pera) [G(2)G(z>) a bee (GG) at oe 


which agrees with (26). 
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The important case represented by Fig. 2, in which 
the Gq channel is continuous, is also amenable to the 
Bode-Shannon approach if a mathematical artifice 
which will be termed a reciprocal sampler is introduced. 
With this device it is permissible to rearrange Fig. 2 
into the form of Fig. 5. 


oe ee. | 
a ° 


ReciPROCAL 
SAMPLER 


> (tT) 


Fig. 5. 


There should be no concern as to the realizability of 
the reciprocal sampler, since the Ga channel has no 
realizability constraints imposed upon it. Apparently 
the transfer function of the reciprocal sampler is equal 
to ¢,,r,(S)/,,-,(2). The optimum realizable filter can 
now be formulated by inspection, treating ¢,,,,(z) and 
$;,r,(Z) as inputs, and considering the equivalent de- 
sired filter to be given by Ga’ =Gad,,r,(5)/rer,(2). The 
optimum filter without regard to realizability is then 

Goont (s) = Srare(2) Ga'(s, Z). 
rr(Z) 
Imposing realizability constraints on (51) according to 
the Bode Shannon approach, it is seen that 


Ree) 2; 1 [2] ) 
et(Z)IL be-(z) 1h 


It will be perceived that the simple manipulation 
used in deriving Fig. 4 cannot be applied to Fig. 5, with 
the consequence that (46) does not appear to follow di- 
rectly from (52). 

The purpose of this discussion has been to investigate 
the applicability of the Bode-Shannon technique to 


(51) 


(52) 
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pulsed filters, thereby supplementing the mathematical 
derivations with the conceptual appeal offered by this 
approach. 

Author’s Comment: The author appreciates the dis- 
cussion by Prof. D. P. Lindorff. Because of the resem- 
blance between transfer functions and pulse-transfer 
functions, correlation functions and correlation § se- 
quences, spectral densities and pulse-spectral densities, 
the extension of the Bode-Shannon approach to the de- 
sign of digital control systems for minimum mean- 
square error. sequence is fairly straightforward. Follow- 
ing the method of Bode and Shannon, the author has 
also checked the result given by (26). In fact, if one is 
versed in the statistical design theory for continuous- 
data control systems, (26) may be written down by 
inspection of the system block digram. The derivation 
of (26) is presented simply to illustrate the design pro- 
cedure for digital control systems based upon the Wie- 
ner-Kolmogoroff theory. 

Difficulties are encountered, however, when the Bode- 
Shannon method is applied to the more important 
problem of designing digital control systems for mini- 
mum mean-square error. This point has also been noted 
by Prof. Lindorff in his discussion. The operation in- 
volved in (52) of Prof. Lindorff’s discussion seems rather 
difficult to perform, since the term in the operating 
brackets is a function of both s and gz. In view of the 
difficulties which arise in treating more important cases, 
the Bode-Shannon approach is not discussed in this 
paper. On the other hand, the design approach pre- 
sented in this paper is quite general, and can be applied 
to other system configurations than that discussed 
above. This technique has been successively applied to 
digital control systems, with noise or disturbance occur- 
ring at other points of the control loop than in the input 
channel of the system. The main purpose of this paper is 
to introduce a simplified approach for the statistical 
design of digital control systems of various configura- 
tions. 
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Determination of Periodic Modes in Relay Servo- 
mechanisms Employing Sampled Data’ 


H. C. TORNG} and W. E. MESERVEf, sENIOR MEMBER, IRE 


Summary—A difference equation approach is presented to de- 
termine the periodic modes in relay servomechanisms employing 
sampled data. In using the method, together with Izawa and Weav- 
er’s method, the investigator can have all the information concerning 
possible periodic modes of a system. 

This new technique is based on expressing periodic sequences in 
terms of orthogonal functions in the discrete sense and balancing 
the harmonics. 


AMPLED-DATA feedback systems containing a 
S relay and a zero-order hold have been extensively 

investigated. The describing function approach for 
these systems was used by Chow! and Russell? to ap- 
proximate the continuous output for sinusoidal inputs. 
Altar and Helstrom,? Kalman,*® Izawa and Weaver,® 
Izawa,’:* and Mullin and Jury’ employed the phase-plane 
approach to study such systems. 

The describing-function approach serves well to ap- 
proximate the possible oscillations due to sinusoidal in- 
puts, but it remains to be seen whether or not the find- 
ings are exhaustive. The phase-plane approach can 
reveal all possible periodic modes. However, this can 
only be achieved through step-by-step computations 
and is generally only applicable to second-order systems. 

Izawa and Weaver’ developed a simple approach to 
determine all possible periods of cyclic response of a 
relay sampled-data feedback system with given relay 
characteristics and sampling period. 


* Received by the PGAC, November 30, 1960; revised manuscript 
received, March 29, 1960. 

+ Cornell University, Ithaca, N. Y. 

1C. K. Chow, “Contactor Servomechanism Employing Sampled- 
Data,” Ph.D. dissertation, Cornell Unive_sity, Ithaca, N. Y.; 1953. 

*C. K. Chow, “Contactor servomechanism employing sampled- 
data,” Trans. AIEE, vol. 73, pt. 2, pp. 51-64; March, 1954. Also see 
F. A. Russell, “Discussion to Chow’s paper,” Trans. AIEE, vol. 73, 
pt. 2, pp. 62-64; March, 1954. 

3W. Altar and C. W. Helstrom, “Phase-plane Representation of 
Sampled-Servomechanisms,” Westinghouse Electric Corp., Pitts- 
burgh, Pa., Res. Rept. No. 60-94410-14-B; September, 1953. 

4R. E. Kalman, “Phase-Plane Analysis of Non-Linear Sampled- 
Data Servomechanisms,” M.S. Thesis, Mass. Inst. Tech., Cambridge, 
Mass.; June, 1954. 

5 R. E. Kalman, “Some aspects of non-linear sampled-data sys- 
tems,” Proc. Symp. on Non-Linear Circuit Analysis, Polytechnic 
Inst. of Brooklyn, Brooklyn, N. Y.; April, 1956. 

6 K. Izawa and L. A. Weaver, “Relay-type feedback control sys- 
tems, with dead time and sampling,” Trans. AITEE, vol. 78, pt. 2; 
May, 1959. 

7 K. Izawa, “On-off control with periodic sensing device,” Trans. 
ASME, vol. 80, pp. 1459-1464; November, 1959. 

8 K. Izawa, discussion following reference 9. 

9F. J. Mullin and E. I. Jury, “A phase-plane approach to relay 
sampled-data feedback systems,” Trans. ATEE (Applications and 
Industry), no. 40, pp. 517-523; January, 1959. 


In this paper the difference equation approach is used, 
together with Izawa and Weaver's findings, to reveal 
all possible periodic modes of relay sampled-data feed- 
back systems of any order in a simplified way. 

The effect of the presence of a dead zone in the relay 
characteristic is studied. The steady-state response of 
the system corresponding to step, ramp, or certain sinus- 
oidal inputs is also discussed. 


DIFFERENCE EQUATION REPRESENTATION 


The system to be considered is shown in Fig. 1. The 
linear portion of the system, with transfer function 
G(s), has p poles. The relay is defined by the following 
relationship 


fr= 1 Cee’ 
= (0 | en| <e 
= — | én <S —e (1) 


r{t) + if 


Fig. 1—A relay sampled-data control system. 

The difference equation relating c, and f,, the system 
output and the relay output of the system at the uth 
sampling instant respectively, can be written as 


Q0Cn+p ote Q1Cn+p—1 ae Q2Cn+p—2 +2 ra Ae i ApCn 
a bo n+p =e Vilatet = bofn+p—2 s= a = AES (2) 


The derivation of (2) can be accomplished either by 
the direct integration or from z-transform tables. The 
order on the right-hand side is at most equal to that of 
the left-hand side. This is due to the presence of a zero- 
order hold in front of the relay. 


PHYSICAL CONSIDERATIONS 


_ When a periodic mode of period NT, where T is the 
sampling period, exists in a relay sampled-data control 
system, the output of the system at sampling instants 
will repeat N values successively; N is an integer. 
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As the periodic output is fed back and sampled, the 
output of the zero-order hold at the sampling instants 
will also repeat a set of N values. Such a set of numbers 
is defined as a sequence. This sequence, when passed 
through the relay, is quantized into a periodic one which 
consists of terms whose magnitudes can only be +1, 0, 
and —1. That is to say, the output of the relay at 
sampling instants can be represented as a periodic se- 
quence, composed of values 1, 0, and —1. 


EXPANSION OF PERIODIC SEQUENCES IN TERMS OF 
ORTHOGONAL FUNCTIONS IN THE DISCRETE SENSE 


In investigating the oscillations that might occur ina 
conventional feedback control system, the possible 
periodic output of the system can be expressed as a 
Fourier series with an infinite number of terms. This 
often makes the complete evaluation impossible. 

In the present study a repetitive discrete sequence, 
representing the output of the system at sampling in- 
stants, is to be determined. One question which will be 
raised is whether or not an expression with a finite num- 
ber of terms can be found to represent the discrete se- 
quence. 

It has been established that this can be done in terms 
of orthogonal functions in the discrete sense. A set of 
functions S,(i) - - - S,(z) with integral argument, 7, de- 
fined when 


Ny <1 <n 


is said to be orthogonal in the discrete sense, if 


Ss SidSa() =0 Lm. 


i=n}1 


(3) 


A periodic sequence P(t) of period N can be written 
in one and only one way as a linear function of 
Si, So, + °°, Sy with constant coefficients.'° 

Trigonometric functions with integral argument form 
an orthogonal set. In terms of this set of orthogonal 
functions, the periodic sequence can be expressed as 
follows: 


k 


ni eee 
20> («. oa + 6, sin n W 


n=0 


if N = 2k + & 


Z 2ri Pete : 
Pj => («. cos 1 ont + 6, sin Oe + ay41 COS Tt 


n=0 


if N = 2k + 2, (4) 


10 T, Fort, “Finite Differences and Difference Equations in the 
Real Domain,” Oxford University Press, New York, N. Y.; 1948. 
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where 

Dh ies 2rm Loe 
Qn = — >. P(m) cosn ii P(m), 

Quel 2am je 
= >> P(m) sin n > O41 =— >, Pim) cos xm, 

N m=0 N N m=0 
nm=1,-:--k, (5) 


DETERMINATION PROCESS 


If an oscillation of period NT prevails in a relay 
sampled-data control system, the output of the system 
will repeat a set of N values. Its output at the mth sam- 
pling instant can be written as 


21n 21n 


k 
Cn = >| v, cosy — + gr siny — |+ t%41008 7m, (6a) 
r=0 N N 


where N =2k+2; and 


27n =| 
Cn = V, COS 7 —— + G, SIN 7 — |], 
dX | N : N 


r=0 


(6b) 


where N=2k-+1. 
The output of the relay at the mth sampling instant, 
fn, can be written as 


k 2rn a 2a7. 
Fn = >| « COs r + h, sin r — | + gr41 COS 7, (7a) 


N N 


r=0 


where N =2k+2; and 


“? 2an _ 2n 
>| « CO P= 5 IS 
N N 


r=0 


where N=2k-+1. 

In determining the possible periodic modes, the value 
of N is first assigned. One question that will be raised 
here is that, since N is arbitrarily assigned, the task of 
identifying the periodic modes will be endless. The fact 
is that the possible values of N for a relay sampled- 
data system is limited and can be readily evaluated by 
Izawa and Weaver’s approach.®® They are determined 
by the half-width of the dead zone and the sampling 
period 7. In an example to follow, the possible values of 
N are only 2, 4 and 6.8 

Once the value of NV is determined, the periodic out- 
put sequence of the relay can be postulated. It seems 
that even with the value of N specified, there are numer- 
ous possibilities for such a postulation. It is actually not 
the case. A systematic approach is presented in the 
Appendix to determine possible sequences that have to 
be considered for a specified NV. 

The postulated sequence, by means of (5), will deter- 
mine the coefficients go, £1, °° * » Zkr Sktty M1 ha, + > + Ie 


in (7a) and (7b): 
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In Appendix I, sets of values of g0,--+, x41, M1, 

- , hy are calculated corresponding to possible se- 

quences. They are ready for universal use as long as the 

investigator normalizes the output of the relay to +1 
and —1. 

Substituting (6a) and (7a) or (6b) and (7b) into (2), 


and collecting terms, one obtains 


2nr 2rr 
s ipa iE Soe amr qr Sere 7 


l=0 


ar 2ar 
= oy Os fs cos 1 + hi, sin a 7 


1=0 


t=) 


. der ea j 
a , cos —— 1 — », sin —— 
eae N N 


l=0 


pe aS aes che (8) 


It should be pointed out here that these N equations 
are linear. Furthermore, for each 7, v, and q, are related 
in two and only two equations. They can be easily 
solved. The determination of the NV unknowns enables 
the periodic output sequence to be evaluated from (6). 
The postulated output of the system at the sampling 
instants will be ¢1, cx, - + - , cy. This is done by assigning 
values of 1 in (6) to be 1, 2, - - - , N respectively. 

The final step is to feed ¢;, co, -- +, cy back to the 
input. If the relay output sequence so produced is the 
same as that upon which the coefficients in (7) are de- 
termined, then there is certainly a periodic oscillation 
which is defined by ¢c1, C2, - + + , wn. 

If the relay output sequence so produced is not the 
same as that postulated, then it is decisively clear that 
there is no such oscillation possible in the servo system. 


ILLUSTRATIVE EXAMPLE 


The system shown in Fig. 1 is to be studied, in which 
the sampling period 7 =1 second, G(s) =1/s(s+1), and 
the relay is ideal. (This is the system studied by Chow,} 
Kalman,’ and Mullin and Jury.®) Even though the 
problem here discussed is of second order, it should 
again be emphasized that the method is equally applica- 
ble to systems of any order. 

The difference equation relating c, and f, can be 
found by direct integration or from a z-transform table 
to be 


Cnz2 — 1.368¢n41 + 0.368c, = 0.368f,41 + 0.264/,. 
That is 
a =1 a; = — 1.368 ad, = 0.368 
bp = 0 b; = 0.368 b. = 0.264, 


p=2 (10) 
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The first task is to determine possible values of JN. 
In the system to be studied, N can only be 2, 4, and 6.° 
D)MieL aN Barc 

From (10), it can be seen that 


p 
Ss ap—-l = 0 


1=0 


It should be remarked here that 


Dp 
‘P OF Sh a 0 
1=0 


is not a coincidence. It is due to the fact that the trans- 
fer function of the plant has a factor 1/s. 

Referring to Table II in the Appendix, the possible 
periodic output sequence of the relay can only be 1, —1. 
The two coefficients gp and gi are 0 and 1 respectively. 

The output of the relay at the mth sampling instant is 


(11) 


As N=2, the output of the system at the nth sam- 
pling instant can be written 


fn = COS TH. 


(12) 


Cn = Vo + 01 COS TH. 


From (10); Di wd) eand (8), the following two 
equations are obtained: 


0 (0.368 — 1.368 + 1)x = 0 
r=1 (0.368 + 1.368 + 1), = 0.264 X 1 — 0.368 X 1 
= — 0.104. (13) 


I 


Y 


Yo and v; are evaluated from (13). Eq. (12) can be 
written as 


Cn = V9 — 0.038 cos an. (14) 
Eq. (13) shows that vp can be chosen arbitrarily, but 
it is not actually the case; v) should be chosen in such a 
way that the output of the system, when fed back to the 
input, will produce the postulated relay output se- 
quence | s=1: 
If vo is so chosen that 


—0.038 < 1) < 0.038, (15) 


the condition to produce the posutlated periodic relay 
output sequence 1, —1 will be met. This shows that the 
system can have the oscillation as defined by (14) 
and. (15). 


2) Let N=4. 


From Table II in the Appendix, the only possible 
periodic sequence of the relay is 1, 1, —1 and —1. The 
four coefficients go, £1, g2 and hy are 0, 1, 0 and 1 respec- 
tively. 


ss 
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The following two expressions can be obtained 


2a 2a 
fn = cos—n-+ sin — n = cos 
4 4 


' TT 
nm + sin 5 Ms (16) 


I 


Cn 


(17) 


T ar 
Vo + v1 COS = n+ v2 cos7mn + gisin — Nn. 
From (10), (16), (17), (8) and (9), a set of four equa- 
tions is obtained: 


(1 — 1.368 + 0.368) = 0 (18a) 
— 0.6320; — 1.36891 = 0.632 (18b) 
1.3680; — 0.632g: = — 0.104 (18c) 

(1 + 1.368 + 0.368). = 0. (18d) 


From (18d) 
aoa a)s 


Solve (18b) and (18c) for v; and qi (v1 and q: are re- 
lated by two and only two equations, 


v1 = — 0.238 
gi = — 0.352. 


The output of the system at the mth sampling in- 
stant is 


Tv T 
Cn = 09 — 0.233 Be n — 0.352 Skee Nn. (19) 


Considering that c, of (19) should produce the postu- 
lated relay output sequence, it is clear that 


—0.238 < 1 < 0.238. 


The results obtained above check with those obtained 
by Kalman. 


DEAD ZONE 


The presence of a dead zone in the relay characteris- 
tic introduces a new output level 0. The periodic mode 
can also be determined by the above procedure, as will 
be shown by the following example. 

In the illustrative example, if the relay characteristic 
is defined as . 


fn =1 Gy, Z2AWAl 
=0 Page g(a 
sae ah én &. — 0.1. (20) 


The possible values of NV can be obtained from the 
Izawa and Weaver method. 
Let N=4; 


Tee et! (postulated). 
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_ Since V =4, there are 4 constants go, g1, g. and /, to be 
determined. They could be determined by using (5) and 
the postulated f,; but, by inspection, it is easily known 
that f, can be written as 


i 
2 = Swe 
2 


(21) 


The output of the system at the mth sampling instant 
can be written as 


Tv Tv 
Cn = Vp + 01 COS a n + v2 cos mn + gi sin mt n. (22) 


From (10), (21), (22), (8) and (9), the following set of 
equations is obtained: 


(1 — 1.368 + 0.368) = 0 (23a) 
— 0.63201 — 1.36891 = 0.368 (23b) 
1.3680; — 0.63291 = 0.264 (23c) 

(1 + 1.368 + 0.368)v2 = 0. (23d) 


From (23d), 


Ve = 0. 
Solve (23b) and (23c) 
v1 = 0.057 
qi = — 0.295, 


Eq. (22) can then be written as 
Tv T 
Cn = to + 0.037 tee e101 295 sine nN. 


The arbitrary constant vp can be determined by con- 
sidering the relay characteristic, shown in (20), and the 
postulated relay output sequence. It is found that when 


—0.043 < vm < 0.043, 


the postulated relay output sequence can be reproduced. 


STEADY-STATE RESPONSE DUE TO VARIOUS INPUTS 


The examples above illustrate the procedure to de- 
termine the periodic modes of relay sampled-data sys- 
tems due to self-excitation. This approach can also be 
applied to systems subjected to step, ramp, or certain 
sinusoidal inputs. This is accomplished by adjusting the 
constant “v.” to ensure that the postulated periodic 
relay output will prevail. 


Step INPUT 


If a step input of magnitude m is applied, the possible 
periodic system output sequence can be expressed as 


piel ( roe e) 
— v, COS —— + Q, Sin — 
C m Vo N qd N 


r=1 


(24) 


+ U,%41 COS TH. 
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Suppose in the system studied in Examples 1 and 2, a 
step input of magnitude 1 is applied, then from (19) and 
(24), the possible steady-state response of the system 
will be 


T T 
Gn = 1 + v9 — 0.238 coher 05352 sit E: nN, 


where —0.238 <u <0.238. 
This result checks with that obtained by Mullin and 
Jury. 


Ramp INPUT 


If a ramp input of “mt” is applied where m is a con- 
stant, the possible periodic system output sequence can 
be expressed as 


feat > ( dan fe =") 
Cn = mn v V, COS — 7, sin ——— 
Ria koe 


r=1 


+ Up41 COS TH. 


(25) 


CERTAIN SINUSOIDAL INPUTS 


When a sinusoidal input with a period P= 7'/b, where 
b is an integer, is applied, the system is essentially ex- 
cited by a step input. If the following relationship exists, 
namely Pb=T7Tw, then the system is excited by a repeti- 
tive sequence of period w. The sequence can be ex- 
panded in the same fashion and the term “gv” can be 
replaced by this expansion and vp. 

The procedure can only be carried out when there is 
an arbitrary constant available. That is to say, the 
transfer function of the linear plant has a factor 1/s. 


CONCLUSION 


A unique difference equation approach is presented 
to determine the periodic modes in a relay sampled- 
data control system of any order. This new technique is 
based on expressing periodic sequences in terms of 
orthogonal functions in the discrete sense and balancing 
the harmonics. It enables the investigator to determine 
decisively whether a specific periodic mode exists or not, 
and if it does, to identify it exactly. In using this 
method together with Izawa and Weaver’s method, the 
investigator can have all the information concerning 
possible periodic modes of a system. 

This approach is extended to the case where the relay 
has a dead zone in its characteristic. 

This approach is also extended to determine the 
possible steady-state response when the system is sub- 
jected to a step, a ramp or a certain sinusoidal input. 

This approach is equally applicable to systems with 
different forms of feedback other than unity. If an ele- 
ment H(s) is in the feedback path, it can be combined 
with G(s) to form G(s) =G(s) H(s). The same procedure 
can then be applied. 
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APPENDIX 


THE DERIVATION OF POSSIBLE PERIODIC OUTPUT 
SEQUENCES OF AN IDEAL RELAY 


Let the number of output levels of a relay element be 
y and the period of the oscillation be NT; the number of 
possible periodic output sequences M of the relay will 
be y%. 

In the case of an ideal relay, y=2, M=2%. When N 
increases, it seems that M will be prohibitively large. 
However, the cases that have to be investigated are 
quite limited, as shown in Table I. 


TAB WE el 
Number of 
Number of sequences to be 
N Number of possible] sequences that considered when 
periodic sequences have to be the plant transfer 


function has an 
integrating factor 


investigated 


Y 4 1 1 
4 16 2 1 
6 64 5 z 
8 256 16 5 


The derivation of Table I and the approach to ex- 
pand it can be presented by considering the case, NV =6. 
M, the number of possible periodic sequences, = 2°= 64. 

These 64 cases can be classified into the following 
groups: Each sequence has 


1) Three (+1)s and three (—1)s 


Quon 4 
= D0 Cases: 
3! 
2) four (+1)s and two (—1)s 
6.5 
=) 1D teases: 
2) 
3) two (+1)s and four (—1)s 
6.5 
— = 15 cases, 
2! 
4) five (+1)s and one (—1)s 
6 
— = 6 cases, 
1 


5) one (+1) and five (—1)s 


6 
———) 0) Gases 
1 ) 


6) six (+1)s 


7) six (—1)s 
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The six cases in group 4) are 


‘le Ps, ee es Bes Ree 
=) eae 


By examining them carefully, one can see that, so far as 
periodic modes are concerned, they actually represent 
one sequence, say 1, 1, 1, 1, 1, —1. It is here defined 
that such a sequence is called a seed sequence. It is 
obvious that for an oscillation of period N, a seed se- 
quence represents NV possible sequences. 

The same argument can be applied to group 5). 
Furthermore, group 5) can be obtained by replacing 
+1 by —1, —1 by 1 in group 4). It is defined here that 
group 5) is a reflection of group 4). Actually groups 4) 
and 5) present just ove case to be investigated. 


Aue, gt ey Ey eee (26) 


Groups 6) and 7) represent actually an oscillation of 
period 1. 

As group 2) is the reflection of group 3), only one has 
to be considered, say group 2). This group presents two 
seed sequences. They are 


(27) 


Groups 2) and 3) present only two cases as shown in 


(27). 
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Finally, let us examine group 1), which consists of 
three seed sequences 


LO DS weaned fa dc raf (28a) 
Teer te eee tie ee (28b) 
2G Eo hy SP eee (28c) 


(28b) is the reflection of (28c) and one seed sequence of 
period 2 


eS ts eet 


It can now be concluded, from (26)—(28), that the 
number of cases to be investigated is five. 

If the plant transfer function has an integrating fac- 
tor 1/s (this is not at all uncommon), then the sum of the 
coefficients on the left-hand side of (2) is zero, that is 


Pp 
>) G1 = 0. (29) 
=) 
Let r=0 in (8); one obtains 
Pp p 
Vo a Ant = Lo Ds (per (30) 
0 1=0 


Considering (29) and (30), one concludes that if 
p 
Ss bp—1 ia 0, 
1=0 


then go has to be zero, that is to say, the sum of the 
terms in the sequence to be investigated has to be zero. 
Under this condition, only two cases from (28) have to 
be considered. 

Table II shows the seed sequences that have to be 
considered for a specified NV. The coefficients of (17) are 
computed by using (5) for each sequence. As soon as the 
investigator normalizes his relay output levels to 1 and 
—1, they can be readily used in (8) and (9). 

Tables I and II are universal in usage and can be 
easily expanded when it is necessary. 


TABLE II 
Coefficients of Expansion 
N Seed Sequences i = A i r 53 5 a 
Agee fol 0 1 
eat (ieee tee 1 0 1 0 1 
pee beh l. e— tl 0 Ps} 0 1/30 2//3 0 
: (LI i seca a ra | 0 1/3 i —1/3 1/3 1//3 
ao ipa pe tee er eet) il 0 2 0 1/2 0 40/2 + 1) 0 z(/2 — 1) 
(Pecos ee ee Ue 1 0 / 2/4 1/2 —4/2/4 1/2 £(4/2 — 2) 1/4 | 0/2 +2) 
£5 91] Me We SII el cet pried ee cea 0 /2/4 1/2 —4/2/4 1/2 4(/2 + 2) 1/4 | 4(,/2 — 2) 
j 1 1 
DB ae iwc Leb mee Ig —1, —1, -1 0 0 1 0 0 a 0 Va 
1, 1 on 14 tit 0 1/2 0 1/2 0 4(4/2 — 1) 0 4(/2+4+ 1) 


Discussion 


A. R. Bergen (University of California, Berkeley). 
The authors present a method for determining whether 
a specified oscillation is possible in a relay sampled- 
data system. A particular periodic relay operating 
schedule is first assumed, and the steady-state sampled 
plant output is then computed; if this output causes the 
relay to operate as assumed, the oscillation is possible. 
Basically, the crux of the problem is to specify the 
steady-state sampled plant output when a _ periodic 
plant input is specified. The authors describe a method 
for accomplishing this; an alternate method, which is 
perhaps more direct, is described below. The method is 
entirely analogous to that used by Bass! for continuous 
relay systems. 

To explain the method consider the authors’ illus- 
trative example, part 1). The output of the relay is 
assumed to alternate between —1 and +1 during suc- 
cessive sampling intervals. Then, assuming the authors’ 
starting phase, the Laplace transform of the relay out- 
put is 


—11-2%147° 1 
is) = 


Ss 1—<;°? Ss 


where z=e*". 

If oscillations in this mode are possible at all, plant 
initial conditions may be found so that the plant output 
assumes its steady-state behavior immediately, without 
a transient. In this case, and assuming a periodic oscil- 
lation, the Laplace transform of the sampled output is 


COG rn Clee Picea ce bee 


Cot cys} Go + c27} 
Prete) 


(2) 


Now if the response to the excitation given by (1) is 
an oscillation in the form of (2), sith cy positive and c; 
negative, the relay will close as assumed and the oscil- 
lation is possible. 

In computing C(z) the initial plant conditions must be 
considered. Returning to the differential equation gov- 
erning the plant 


eté=f. (3) 
The Laplace transform is found to be 


s(s + 1)C(s) = F(s) + (s + 1)co + Oo. (4) 


1R. W. Bass, “Equivalent linearization; nonlinear circuit syn- 
thesis and the stabilization and optimization of control systems,” 
Proc. Symp. on Nonlinear Circuit Analysis, Interscience Publishers, 
Inc., New York, N. Y., vol. VI, pp. 163-198; April, 1956. 
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Inserting F(s) from (1), 


—1 1 — <2} Co Co s 
C(s) = = aaa 0) 

s(s +1) 1+27 s  s(s+1) 

Using a table of z transforms, 
—Tz-} (1 as e-T)g-1 
Ce) = - = =e 
Q—-s)a4+s45 Ud+s)0—e%27) 
: 1 235 (SER peel 

| Co Col Cale (6) 


1—;, (A —2z (1 — e774 ' 


The question is now this: can (6) be put into the form 
of (2), with co positive and c; negative? If so, the oscilla- 
tion is possible. Equating (2) and (6), 
et Ca ee) 


Cot Oe) = 


if et gl 
pe EES Pay sae ee 


1 —@¢ Tz! 


a condition which is clearly impossible, unless the de- 
nominator of the third term on the right of (7) is can- 
celled. This occurs when 


1—et 


iat toes 


Co 


; (8) 
and in this case, 


1—er 
cotcagzt=cot+ (« —T+2 os) peeled | AE: 


The authors assume 7'=1. Then 


Co + c1z-! = Co + (co — 0.076)27}. (10) 

Now if 0<co<0.076, co will be positive, c; will be 
negative, and the postulated oscillation is possible. This — 
same condition is obtained by the authors. 

The applicability of the method is obviously not re- 

stricted to the illustrative example. The general testing 
testing procedure remains the same. Some of the arbi- 
trary initial conditions are used up in satisfying the re- 
quirement of the assumed periodicity. The remaining 
initial conditions are used in attempting to satisfy the 
requirements of sign of the c;. 
f The method of the authors and that of this discus- 
sion share the advantages of being exact and of not being 
limited by dimensionality. Both methods, however, 
require an exhaustive investigation to reveal possible 
modes, and furnish necessary but not sufficient condi- 
tions for their existence. 

Authors’ Comments. The authors believe that the 
approach we present furnishes a simple and exhaustive 
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study of the periodic modes of a relay sampled-data 
system. The relay output schedule is not entirely as- 
sumed, but is based on the work by Izawa and Weaver. 
It is the authors’ belief that this method, together with 
the Izawa and Weaver approach, will be a powerful 
tool in investigating the periodic modes of a relay sam- 
pled-data system. 

If one cares to look carefully at the development of the 
art of sampled-data control systems, he will realize that 
it has partly been an effort to bring to the discrete do- 
main the techniques useful in studying the continuous 
control systems. In linear system analysis, for example, 
the prevailing transformation methods, and certain 
processing techniques, all represent attempts in this 
direction. Correspondingly, in nonlinear system analy- 
sis, we have the phase-plane approach and describing 
function approach. But in doing this, one has to be 
aware of the fundamental differences between a con- 
tinuous system and a discrete system. For instance, 
certain oscillations in a sampled-data relay system will 
not be found in a continuous system with the same plant 
and relay characteristics. The remark, “The method is 
entirely analogous to that used by Bass! for continuous 
relay systems,” would be quite justified by casually 
glancing through the opening paragraphs of Bass’s 
work, and, in the phrase, “The discovery of all periodic 
solutions of a given system of ordinary non-linear dif- 
ferential equations,” substituting the word “difference” 
for “differential.” But one can not ignore the basic dif- 
ference between the nonlinear differential equation and 
a nonlinear difference equation. If one does examine in 
detail Bass’s work and other more recent work? in 
Russia on the continuous relay system, he will realize the 
significant difference in analyzing a continuous system 
and a discrete system. It is the authors’ opinion that 
Bergen’s remark is irrelevent. 

We find that Bergen’s alternate method is in essence 
really “entirely analogous” to our method. We are fully 
aware of the usefulness of the z-transform technique, 
but we do not believe that its application in the present 
study will make it more direct, as asserted by the dis- 
cussor. 

Let us now look at this alternate method. 

First of all, our starting phase in the illustrative ex- 
ample is +1, —1 instead of what Bergen “assumes” it 
to be. 


2M. A. Aizerman and F. R. Gantmakher, “Determination of 
periodic modes in systems with piecewise-linear characteristics com- 
posed of segments parallel to two specified straight lines,” Automation 
and Remote Control, vol. 18, pp. 105-120; February, 1957. 
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For convenience of comparison, we adopt the —1, 1 
starting phase in the following discussion. 


Bergen’s Method 
1) The relay output 


Y(s) = -— as — . 
s 1+ 21 

One has to take the Laplace 
transform of a discontinuous 
function, and this needs ma- 
nipulation to obtain the closed 
form. When JN isa large num- 
ber, this task is somewhat 
formidable. 


2 


wa 


The system output at the 
sampling instant: 


(2) =e toast: 
cot az! 
(Ui a ire) 


Manipulation is needed to 
obtain the closed form. 
3) The plant equation: 

Write down the differential 
equation from the transfer 
function of the given plant. 
Take the Laplace transform. 
This step is to bring in one 
more unknown, ¢o in this case, 
more in higher order systems. 


s(s + 1)(s) 
= V(s) + (s+ leo + Co. 


4) Take the z-transform and ob- 
tain a bulky equation for 
C(z). This could be more 
cumbersome with a large 
period of oscillation. 

Equate the expression in step 
2) to the long expression in 
step 6) and arrange terms. 


5) 


wm 


6 


we 


Inspect the equation in step 
5), attempt to find suitable 
condition for éo. It should be 
remarked here for a higher or- 
der system there seems to be 
no general rule to follow to 
eliminate these unnecessarily 
introduced unknowns. 

7) Find the range for co and de- 
duce c by keeping the valid- 
ity of the result in mind. 


Our Method 


Y(n) = — Cos ur 


nm from O0—> ©. 


This can be obtained by inspec- 
tion in simple cases, or computa- 
tion from (5), or referring to 
Table II in the Appendix. 


c(n) = vo + 21 Cos rn. 


No manipulation is needed, just 
write down the expression from 


(6). 


From the transfer function of the 
given plant, the difference equa- 
tion is obtained by direct integra- 
tion or with the help of a g-trans- 
form table. 


Gaia = 1 308c gt 0. 308e 
= 0.368yn41 + 0.264yn. 


No additional unknown is intro- 
duced. 
No corresponding step necessary. 


Substitute coefficients obtained 
from steps 1), 2), and 4) to 8) to 
obtain two linear equations in 
terms of vp and % which can be 
easily solved. 

No corresponding step necessary. 


Check the validity by merely 
inspecting the signs of co and 4. 


After the above step by step comparison, Bergen’s 
claim that his method is perhaps more direct clearly 


needs qualification. 


We are grateful to Professor Bergen for his discussion. 
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Analysis of Pulse Duration Sampled-Data 
Systems with Linear Elements 


R. E. ANDEENT, MEMBER, IRE 


Summary—This paper deals with the application of z-transform 
techniques to the analysis of sampled-data systems in which signals 
appear in pulse duration modulated form. The characteristics of 
pulse duration sampled data are first briefly described. An analyt- 
ical means for studying the transient response and stability of sys- 
tems which use this kind of data is then presented. Finally, a com- 
parison is shown of analytical results with tests on several repre- 
sentative experimental systems. A short discussion of the relative 
advantages and disadvantages of pulse duration modulation in con- 
trol systems is also included. 


INTRODUCTION 
ap HE STUDIES of sampled-data systems made up 


to now have ordinarily assumed that sampling is 
a kind of pulse amplitude modulation process. The 
output of the sampler is usually thought of as a sequence 
of equally spaced, equal duration pulses whose ampli- 
tudes are proportional to the input to the sampler at the 
sampling instants (constant sampling frequency and 
constant pulse width). There are, however, other ways 
of representing sampled intelligence by a modulation 
process. One of these requires that the sampler produce 
an output pulse train made up of equally spaced, con- 
stant amplitude pulses, whose individual pulse dura- 
tions vary as a function of the input to the sampler at 
the sampling instants. In contrast to the former, this 
is a pulse duration modulation process. These two types 
of sampled data illustrated in Fig. 1 will be referred to in 
this paper by the terms pulse amplitude sampled data 
(PASD) and pulse duration sampled data (PDSD).1 
The classical example of a system which employs 
PDSD illustrated in Fig. 2, was described by Gouy 
approximately sixty years ago.? This is a temperature 
regulating system which consists of an oven heated by 
a resistive heating element, a source of electrical poten- 
tial, a relay, a mercury thermometer, and an electrical 
contact mounted on a motor driven eccentric. When- 
ever the electrical contact becomes immersed in the 
mercury, current flows in the coil of the relay which 
interrupts the heating current. Since the contact is 
periodically dipped into the mercury, the heating cur- 
rent consists of a sequence of pulses. These pulses have 
equal amplitudes, but their widths depend upon the 
time that the contact is immersed in the mercury. The 


* Received by the PGAC, August 13, 1959; revised manuscript 
received, April 18, 1960. 
a t Sperry Phoenix Co., a division of Sperry Rand Corp., Phoenix, 

abe, 

*H. S. Black, “Modulation Theory,” D. Van Nostrand Co., Inc., 
New York, N. Y.; 1953. 

* M. Gouy, “On a constant temperature oven,” J. Physique, vol. 
6, ser. 3, pp. 479-483; 1897. 
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Fig. 1—Comparison of sampled data. (a) Input; 
(b) PASD output; (c) PDSD output. 


time that the contact is immersed in the mercury de- 
pends upon the level of the mercury in the thermometer, 
and the level of the mercury depends upon the tem- 
perature of the oven. Thus the current pulse widths are 
varied in proportion to the temperature of the oven, and 
feedback control of the oven temperature is attained. 
One of the reasons for the past emphasis on PASD is 
that the mathematical methods for analyzing control 
systems based on this type of data have been extensively 
developed. A pulse amplitude sampler is essentially a 
linear device in the sense that inputs and outputs are 
governed by the superposition principle. This is not the 
case with a pulse duration sampler, and as a consequence 
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LI 


Fig. 2—Gouy temperature regulator. 4 =motor driven eccentric, 
B=electrical contact, C=relay coil, D=normally closed con- 
tacts, E=mercury thermometer, F=resistive heating element, 
G=oven. 


the development of mathematical means for analyzing 
PDSD systems has been hindered. On the other hand, 
PDSD systems have a number of advantages which 
recommend their use. This paper shows how the 
Principle of Equivalent Areas? may be applied to the 
mathematical analysis of PDSD systems with linear 
continuous elements. 


THE CHARACTERISTICS OF PULSE DURATION 
SAMPLED DATA 


The PDSD described in this paper is generalized 
somewhat to include the possibility of a negative sam- 
pler output as well as positive and zero outputs. The 
polarity of the sampler output depends upon the po- 
larity of the input to the sampler. Furthermore, it is 
assumed that the sampler output pulses all have their 
leading edges at the sampling instants t=nT. The loca- 
tions of the trailing edges of the pulses are assumed to be 
a linear function of the input to the sampler at the 
sampling instants. In the original Gouy temperature 
regulator the output pulses are symmetrical with respect 
to the sampling instants, and consequently differ 
slightly from the assumptions which we have made for 
mathematical convenience. However, since the tem- 


3R. E. Andeen, “The principle of equivalent areas.” (To be pub- 
lished in Trans. AIEE.) 
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perature change is assumed small during a sampling pe- 
riod this difference is of little consequence. 

Consider Fig. 1. If the input to a pulse duration 
sampler is r(t), then the sampled output rp(é) is given 
by the following expression: 


rp(t) 
=) fi = sign r(mT)[u_ilt — nT) — ult — nT — hn)|, (1) 


n=0 


where h,, the variable pulse duration, is equal to the 
following 


ln = a| r(nT) | 
h, = T 


a| r(nT) | Sey 
a| r(nT) | il: 


It is apparent at the very outset that while a pulse 
amplitude sampler is a linear element in the sense that 
inputs and outputs are governed by the superposition 
principle, a pulse duration sampler is a nonlinear ele- 
ment. Two types of nonlinearity are present. First, 
superposition does not apply; that is, the sampled sum 
of two inputs does not equal the sum of two sampled 
outputs corresponding to these inputs. Second, the 
sampler output saturates when a| r(nT) | = T since obvi- 
ously the pulse duration cannot exceed the sampling 
period. In the following analysis we will always assume 
that the input to the sampler is not greater than the 
value which causes saturation. In other words, we will 
restrict the sampler operation to a quasi-linear un- 
saturated region. 

Systems with PASD usually depend upon the use of a 
hold circuit to convert the sampler output pulse train 
into a staircase type signal. This is necessary for two 
reasons. First, converting the pulsed input to a stair- 
case input allows more energy to be delivered to the 
system during a sampling period. If one increases the 
pulse amplitude instead, one runs into the problem of 
saturation. Second, this reduces the amount of ripple in 
the output because the hold circuit acts as a low pass 
filter on the pulsed signal. Systems with PDSD must 
depend on the system elements themselves to perform 
the filtering function. These systems must include some 
energy storage device such as a motor inertia or an oven 
heat capacity which will smooth out the effects of the 
pulsed nature of the PDSD input. To provide adequate 
filtering, the following discussion and analysis conse- 
quently assumes that a pulse duration sampler always 
works into a combination of system elements which in- 
cludes a low pass element having a largest time constant 
at least twice the sampling period (or if the roots are 
imaginary, the reciprocal of the real part of the largest 
root is at least twice as great as the sampling period). 

There is a great deal of similarity between the output 
of a pulse duration sampler and the output of an on-off 
relay controller. Both outputs are pulsed and may 
assume only the amplitudes +E, zero, and —E. The 
principal difference between the outputs is that the 
sampling period of the pulse duration sampler is fixed, 
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while both the pulse width and pulse repetition rate of 
the on-off relay controller depend upon the input to the 
controller. This similarity is not exploited in the fol- 
lowing. 

The following definitions will be important in the 
analysis of PDSD systems. We assume that the absolute 
value of the input to the pulse duration sampler, 
| (2) | , is bounded and Rnax is equal to the largest value 
of | (2) | which we measure. The maximum pulse width 
aRmax by assumption is less than, or equal to, the sam- 
pling period 7. The modulation ratio, 8, is defined as 
the ratio of @Rmax to the sampling period 7’, 


ORmax 


ee a (2) 


B must lie between zero and one. The filter ratio, a, is 
defined as the ratio of the sampling period to the largest 
time constant of the system, 


E 
dir ax 


a= 


(3) 
By assumption, a@ must lie between zero and one half. 


ANALYSIS OF LINEAR SYSTEMS 


This analysis is based on the principle of equivalent 
areas (Appendix) and the z-transform method. The 
principle of equivalent areas expresses the conditions 
under which two input signals to a dynamic element 
produce essentially the same output. According to this 
principle, two input signals are dynamically equivalent 
if their integrals evaluated over corresponding sampling 
intervals are equal, and the sampling frequency deter- 
mined by this interval is at least twice the highest sig- 
nificant frequency in the frequency spectrum of the 
dynamic element. Expressed mathematically, two input 
signals are equivalent for suitably small T if 


nT nT 
i= f r(t)dt ai 
(n—1)T (n—1)T 


for all values of 7 sufficient to cover the time interval of 
interest. 


r (fat T 


PULSE 
DURATION 
SAMPLER 


CONTINUOUS 
ELEMENT 
G(s) 


Fig. 3—System with pulse duration sampler. 


Consider Fig. 3 and assume that r(¢) is some input 
control signal. r(t) is sampled by the pulse duration 
sampler. The output rp(t) is PDSD and is the input to 
the continuous linear element G(s). We are concerned 
with determining the behavior of the output c(#) when 
r(#) and the system parameters are known. To do this 
the pulse duration sampler is replaced by an equivalent 
circuit comprised of a pulse amplitude sampler in cas- 
cade with a hold circuit described by a transfer function 
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H(s). Restricting attention to the sampling instants, 
the z-transform method is applied to the linearized sys- 
tem to determine the equivalent z-transform of the over- 
all PDSD system. Transient response and _ stability 
boundaries can then be investigated by z-transform 
techniques in the usual way. 

The replacement of the pulse duration sampler by an 
equivalent circuit is accomplished as follows. Refer to 
Fig. 4. This figure shows some arbitrary sampler input 
r(t) and the corresponding PDSD output rp(t). A 
PASD signal rp/(£) is constructed from rp(t) according 
to the Principle of Equivalent Areas. 7 is chosen as an 
interval over which area must be conserved. The areas 
of corresponding pulses of 7p(¢) and 7p’(¢) are required 
to be identical. Under these conditions the equivalent 
PASD is related to r(t) by the following expression, 


rol) = —— Dra) 


max n=0 


-[u_ili — nT) — u_ilt — nT — BT)]. (4) 


(a) 


(c) 


Fig. 4—Equivalent sampled data. (a) Input signai; 
(b) PDSD signal; (c) equivalent PASD signal 
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The Laplace transform of (4) is 


me ID 


D> r(nT)e-”?, (5) 


max n=0 


Ra'(s) = 


Eq. (5) can be looked on as representing an amplitude 
sampler in cascade with a hold circuit with the charac- 
teristic 


H(s) ie 
5 igs: ( ) 


—— H(s) G(s) 


Fig. 5—Equivalent circuit for system with pulse duration sampler. 


Fig. 5 thus represents a linearized equivalent circuit 
corresponding to Fig. 4 which includes an amplitude 
sampler in cascade with a hold circuit to represent the 
pulse duration sampler. One must remember that the 
validity of this equivalent circuit depends upon G(s). 
G(s) must act as a low pass filter whose dominant time 
constant is at least twice the sampling period (0<a<$). 
This, as pointed out previously, is usually the case in a 
practical system where filtering action is necessary to 
attenuate ripple in the output. 

The over-all equivalent z-transform c*/r* which 
represents the system comprised of a pulse duration 
sampler in cascade with G(s) is easily determined for any 
particular G(s). Consider the following example. As- 
sume that G(s) is a simple lag network whose transfer 
function is 


G(s) = 


Tare (7) 


The relations in Fig. 5 are represented mathematically 
by 


C(s) = H(s)G(s)R*(s) (8) 
C*(2) = Z[H(s)G(s) |R*(2). (9) 


The z-transform of H(s)G(s) is evaluated by complex 
integration as shown by Jury.* 


1 pote 1 
Z[H(s)G(s)] = sie H(p)G(p) Pagar (10) 
E(1 peice A 
ZH = NaS. pA 
dp 
oe a}. (11) 


4.1. Jury, “Analysis and synthesis of sampled data control sys- 
tems,” roe: AIEE, vol. 73, pp. 332-346; September, 1954. 
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Evaluating by residues we obtain, 


E {| 1 1 | 
Ree Neleeean) [eee a 


1 e ATU-8B) 
— gal — 2 
" iF —gi{t— —a\|t poe 


The following equation results from combining terms 
and substituting a=AT: 


Z| H(s)G(s) 


Z| H(s)G(s)| = 


shes 1 sel Piptas is \ tl 8) 
Ie Iie Sea 1 eae 
Eq. (13) can be reduced to the following, 
Z[H(s)G(s)] = AC ene! Sens, 
bes en? et Ps Rate 
Thus the result is 
C*(z) bees e-*(e%® — 1)g-4 (15) 


R*(z) F Wetopes WW ey 


The results of similar derivations for some PDSD sys- 
tems with different continuous elements are shown in 
Table I. The reader can readily extend this table to 
include any G(s) of special interest by following the 
procedure described above. 


TABLE I 
List oF EQUIVALENT PDSD SystTEM z-TRANSFORMS 


G(s) Equivalent PDSD System z-Transform C*/R* 
A E (So - = 
s+A ener Ma ea 


AB E [( a ) (ev® — 1)e-ve1 
(eit A is B) Reva WaiNOo 17, f =senzns 
vy \ (e% — ey 
0 e _ -) 1—e 41 
Te@(1 — e)g-1 
(1 — e227) 


A E 
Rmax 


Bz? 


1-7} 


ADVANTAGES AND DISADVANTAGES OF USING PULSE 
DURATION SAMPLED DATA 


Some of the advantages of using PDSD in a feedback 
control system are the following: 


1) PDSD is less susceptible to the influence of noise 
which manifests itself by random amplitude vari- 
ations since the information contained by an indi- 
vidual pulse resides in the length of the pulse. Ran- 
dom amplitude variations can be eliminated by 


clipping. 
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2) The output stage of the pulse duration sampler can 
be a simple relay circuit. Thus very high power 
gains can be realized with little equipment com- 
plexity and relatively low cost. 

3) In the presence of sticking friction or backlash it 
is more effective to apply a large torque for a short 
time than a smaller torque for a relatively longer 
time. A system based on the use of PDSD to drive 
an output motor against a load with friction and 
backlash applies full torque to the load for a 
greater or shorter time depending upon the input 
to the sampler. This is consequently more effective 
than a system based on the use of PASD (with a 
hold circuit) which applies a torque proportional 
to the input to the sampler for an entire sampling 
period. In order to obtain an acceptable torque 
gradient (pound-foot/volt of error) with the sys- 
tem based on PASD (or for that matter, also sys- 
tems using continuous data) an output motor is 
often used which has a much higher rated torque 
than that which is actually required for maximum 
load. This is inefficient, more costly, and intro- 
duces the problem of torque limiting. 

4) Torque limiting is particularly important in air- 
craft applications where exceeding the structural 
limits of the aircraft would be disastrous. When 
PDSD is used to drive the output motor, torque 
limiting is easily accomplished since only one 
value of torque is ever applied. The associated 
circuitry can easily be designed to be fail safe. 

5) Compared to an on-off regulator, the pulse dura- 
tion sampling regulator has the following advan- 
tage. With the on-off regulator the steady-state 
ripple in the output is a function of the gain of the 
system and the dead band of the relays. The fre- 
quency of the ripple may change with system gain, 
but the amplitude is fixed by the dead band. With 
the pulse duration sampling regulator, the fre- 
quency of the ripple is a constant which is com- 
pletely under the control of the designer. The 
amplitude of the ripple is still a function of the 
gain of the system, however, by increasing the 
sampling frequency and reducing the system gain 
at the same time to maintain a constant ratio of 
on time to total time, the ripple can be reduced to 
any desired level. 

6) No hold circuit is required. 


Some of the disadvantages of using PDSD are the 
following: 


1) PDSD can only be used where low pass filter ele- 
ments provide adequate filtering to reduce the 
ripple due to the pulsed input to an acceptable 
level. 

2) A pulsed input and consequent pulsating torque 
may cause greater heating and mechanical stresses 
in an output motor and associated load elements. 
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This can be reduced somewhat by introducing 
additional filtering. 


EXPERIMENTAL RESULTS 


A series of tests were made to obtain an experimental 
confirmation of the theoretical results given in Table I. 
A pulse duration sampler and various linear systems 
were constructed with electronic and relay elements 
and an analog computer as described in the Appendix, 
and were subjected to step and sine-squared pulse 


es PULSE CONTINUOUS | 
r DURATION ELEMENT 
2 SAMPLER | ', As) | 


Fig. 6—PDSD system with feedback. 


transients. Figs. 3 and 6 illustrate the configurations of 
the open-loop and feedback systems tested. Several con- 
tinuous elements were considered: 


System I 
G(s) A 
ee a eee 
(Se A) 
System IT 
AB 
G(s) = 
(SAG 4-38) 
System ITI 
A 
~ 35(s + A) 
Three combinaations of parameters were studied, 
Case a) a=4 
B=3 
T=8 seconds 
Case b) a=+ 
p=3 
T =8 seconds 
Case c) a=} 
B=} 
T=8 seconds. 


The theoretical responses which would be obtained 
from the open-loop and feedback systems illustrated in 
Figs. 3 and 6 for the given transient inputs were also 
calculated by means of the equations in Table I and 
the usual z-transform methods. The experimental and 
theoretical results for Case b) are compared in Figs. 7 
and 8 for the open-loop and feedback cases respectively. 
The results for Cases a) and c), not shown here, were 
similar. From these figures it appears that in both the 
open-loop and feedback situations the equations of 
Table I were able to predict the performance of the 
PDSD systems studied to a degree of approximation 
satisfactory for practical purposes. 
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The stability boundaries of the PDSD feedback sys- 
tems illustrated by Fig. 6 were also calculated from the 
equations of Table I and compared with experimentally 
determined stability boundaries. The results are sum- 
marized in Table II. It is apparent that the agreement 
between the experimental and calculated results is good. 


TABLE II 
GAIN AT STABILITY BOUNDARY 
a Ke Ks 
SSS Case Calculated Experimental 

a) 9.34 9.00 
I b) 11.1 9.95 

c) Seal 29.3 
a) 10.47 9.25 

II b) 16.1 1520 

Cc) 41.0 33.8 
a) 1.19 Le20 
Ill b) 0.474 ONS 
Cc) 1.045 1.10 

CONCLUSIONS 


An approximate method of analyzing PDSD systems 
with linear elements can be based on the z-transform 
method if the PDSD signal is replaced by an equivalent 
PASD signal. This replacement is based on the prin- 
ciple of equivalent areas which requires that the areas of 
the corresponding pulses of the actual and equivalent 
sampled signals to be equal. This type of linearized 
analysis is based on two assumptions. First, the opera- 
tion of the pulse duration sampler is restricted to a 
quasi-linear region where input variations can be trans- 
lated into pulse duration variations. Second, it is as- 
sumed that the continuous elements following the pulse 
duration sampler include a low pass element which has 
a time constant which is at least twice as great as the 
sampling period. This condition is not very restrictive, 
for any practical system requires filtering to reduce the 
effects of the pulsed signal. 

When we replace the PDSD signal by an equivalent 
PASD signal we find that a pulse duration sampler can 
be represented by an equivalent circuit consisting of an 
ideal amplitude sampler and a special hold circuit whose 
parameters are a function of the modulation ratio and 
the maximum value of the input. Based on this equiva- 
lent circuit the system z-transforms of a pulse duration 
sampler in combination with several different contin- 
uous elements were derived and are shown in Table I. 
A series of tests which were designed to obtain an ex- 
perimental confirmation of these equations gave satis- 
factory results. We conclude that the method of study- 
ing the behavior of PDSD systems described here, sub- 
ject to the above assumptions, will give results satis- 
factory for practical engineering calculations. 

The success of this approximate method of analysis 
suggests that digital or discrete compensating elements 
could be employed to improve the transient response or 
other characteristics of PDSD systems. A large body of 
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literature already exists which treats the application of 
digital or discrete compensation to PASD systems. 
Much of this can probably be carried over to PDSD 
systems. Mixed PASD and PDSD systems seem fea- 
sible. 


APPENDIX I 
PRINCIPLE OF EQUIVALENT AREAS 


The principle of equivalent areas* is an approximation 
expressing the conditions under which two inputs to a 
dynamic element produce essentially the same output. 
The derivation starts with the convolution integral 
relating the output c(t), input 7(¢), and the impulsive 
response of the dynamic element A(t), 


c(t) = h(r — t)r(t)dt. 

We replace h(t) by a staircase approximation h,(t) on 
the basis that /,(t) changes value at a sampling interval 
short enough so that the frequency spectra H(jw) and 
H,(jw) are similar. A sampling frequency at least twice 
the highest significant frequency of H(jw) will ordi- 
narily accomplish this. Assume that 7 is an integral 
number of sampling periods, 


c= md NM = QL 230 tered, 
Then we can write, 
mT 
c(mT) = hp(mT — t)r(t)dt. 
Furthermore the integral from — 2» to mT can be 


broken up into-a sum of integrals over the intervals 
(n—1)T to nT where n goes from — © to m, 


c(mT) = Lf mmr — into]. 


However, during the interval (n—1)T to nT, h,(mT —2) 
is constant and equal to h,(m7T—nT). Consequently we 


obtain 
r= > a hy(mT — nT)r(t)dt 
c(mT) 22 eee (m Yr (2) i 
c(mT) = 5s h,(mT — nT) ss r(t)dt. 


n=—co 


Now define I(nT) by 


nT 
I(nT) = i r(t)dt. 
(n—1) T 


(m—1)T 


Then we obtain our final result, 


tn) oO” EAT RTE 


n=—0 


APPENDIX II 
EXPERIMENTAL EQUIPMENT 


The PDSD systems represented by Figs. 3 and 6 


1960 


were set up in the laboratory using the following equip- 
ment: 


1) low frequency square wave generator, 

2) electronic analog computer, 

3) 4-channel recorder, 

4) 300-volt de power supplies, 

5) (amplitude) sampler and zero-order hold unit, 
6) pulse duration relay unit. 


The first four items are general laboratory type equip- 
ment, and the last two items were specially built for 
sampled-data studies. 


Sampler and Zero-Order Hold Unit 


This device operates to convert a continuous input 
voltage into a staircase output voltage which is equiva- 
lent to the output which one would obtain by the proc- 
ess of sampling and holding. The basic performance 
characteristics of this device are listed below: 


range of sampling frequencies—0.20 to 2000 cps, 

attenuation—approximately 0.84, 

maximum voltage of input signal— +25 volts peak, 

phase shift—less than 10 degrees at highest frequency 
of the frequency selector setting. 


The operation of the S-H unit is based on the principle 
shown in Fig. 9, but the switching is done electronically. 
The switches shown in Fig. 9 operate once each sampling 
interval. While one capacitor is charging up to the am- 
plitude of the input, the other capacitor is disconnected 
from the input and-supplies the output voltage to the 
load. A buffer amplifier is required to prevent discharg- 
ing the capacitor connected to the load. A complete de- 
scription of the electronic circuit, and details of opera- 
tion and calibration are given by the author.’ 


Pulse Duration Relay Unit 


The relay unit consists of five channels, each including 
an electronic amplifier which operates a sensitive relay. 
When a positive voltage is applied to the channel input 
the relay operates; for a negative voltage the relay drops 
out. Because of the high amplifier gain the differential 
voltage between pull-in and drop-out is insignificant. 


PDSD System Operation 


In setting up the systems described by Figs. 3 and 6 
the continuous elements were simulated by the ana- 
log computer. A pulse duration sampler was con- 
structed from the S-H unit, the relay unit, and an in- 
tegrator and an amplifier of the analog computer. The 
operation of this simulated pulse duration sampler was 
as follows. The continuous input to the pulse duration 
sampler was applied to the S-H unit which converted 
this continuous voltage to a staircase type voltage. The 
output voltage of the S-H unit was applied to the input 


BR E. Andeen, “A wide-range sampler with zero-order hold,” 
Control Engrg., vol. 6, pp. 149-151; May., 1959. 
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Fig. 10—Simplified block diagram of simulated 
pulse duration sampler. 


of an analog computer amplifier, which is designated 
Al in Fig. 10. The input to the integrator was a con- 
stant voltage, and its output was also connected to the 
input of amplifier 41, but with a polarity opposite to 
that of the S-H unit output voltage. At the beginning 
of each sampling period the integrator was reset to zero 
output voltage, so that its output appeared as a saw- 
toothed voltage. A relay operated by the output voltage 
of amplifier 41 completed the pulse duration sampler. 
When the relay operated it applied the constant voltage 
+E to the continuous element which followed it in the 
system. At the beginning of each sampling period the 
relay would first operate, applying +£ to the contin- 
uous element, because the polarity of the S-H unit out- 
put was the proper polarity to cause operation and the 
integrator output was zero. During the sampling in- 
terval the output of the integrator would increase and 
when it exceeded the S-H unit output the relay would 
drop out. It is clear that the time interval required for 
the integrator output to build up to the value of the 
S-H unit output was directly proportional to the mag- 
nitude of the output of the S-H unit. Thus the duration 
of the +£ output from the relay was made proportional 
to the input to the simulated pulse duration sampler at 
the sampling instants. Additional switching was ar- 
ranged to take care of changing the polarity of the out- 
put and the input voltage to the integrator when the 
continuous input to the simulated pulse duration 
sampler changed polarity. Fig. 10 is simplified by show- 
ing the connections for operation with only one polarity 
of the input r(f). 
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The Analysis of Cross-Coupling Effects on the 
Stability of Two-Dimensional, Orthogonal, 
Feedback Control Systems” 


DAVID BRUCE NEWMANT, SENIOR MEMBER, IRE 


Summary—The analysis of two-dimensional cross-coupled feed- 
back control systems is studied because many such systems are in 
existence, and methods for their analysis and criteria for measuring 
their performance have had only limited attention in the technical 
literature [{1|-[6]. Previous literature has generally been in the area 
of matrix analysis of multidimensional or multiloop systems, and 
has not treated the handling of poles and zeros in the complex s- 
plane or cross-coupling, as they are treated in this paper. 

An algebraic model is developed for a generalized system, using 
the assumption that all components of the system, including the cross- 
coupling transfer function, may be represented by expressions in the 
complex frequency domain which have real-time domain equivalents 
describable by linear differential equations. 

The stability of the two-dimensional, cross-coupled, feedback 
control system is studied by complex-plane methods, and a method 
is developed for determining whether a system is stable or not. A 
technique is also indicated to show how an unstable system may be 
made stable by the introduction of appropriate cross-coupling. 


I. INTRODUCTION 


ANY two-dimensional, orthogonal, feedback 
M control systems are in existence today. An ex- 

cellent example of such a system is an auto- 
matically tracking radar system in which the antenna 
is continually directed in each of two orthogonal planes 
of action so as to minimize an error signal. The purpose 
of this system is to keep the antenna pointed towards 
the target, and thus establish the vector to the target 
from the point of observation. 

In two-dimensional systems, cross-coupling occurs 
between the two orthogonal systems in many instances, 
because these systems share a common facility, such as 
a power supply or radio transmission link. The normal 
procedure for the analysis of two-dimensional systems 
has been to assume complete independence of the two 
systems, and to analyze each system independently. 
This paper discusses analytically the effect of cross- 
coupling on the stability of two-dimensional systems. It 
also shows that the introduction of cross-coupling into 
an unstable system can make the system stable. 


II. ALGEBRAIC MODEL FoR ANALYSIS 


To establish a general framework for the analysis of a 
feedback control system it is necessary to develop an 
algebraic model. 


* Received by the PGAC, July 9, 1959; revised manuscript 
received April 19, 1960. 
{ Fairchild Astrionics Div., Wyandanch, Long Island, N. Y. 


A one-dimensional feedback control system is gen- 
erally represented by the following diagram, 


where s is the complex Laplace operator, R(s) is the 
transform of the reference input function and C(s) the 
controlled output function. C(s) is to follow R(s) as 
nearly as possible. R(s) and C(s) are compared contin- 
uously. /(s) is the error function obtained by subtract- 
ing C(s) from R(s), and C(s) is-determined by E(s) 
operating through a physical device whose transfer 
function is G(s). An implicit assumption is that each 
term involved may be represented in real time by a 
linear differential equation, which is transformable from 
the time domain to the complex frequency (s) domain 
by the Laplace transformation, 


FQ) = Ll = [pone (1) 


The total transfer function for the general one-dimen- 
sional system is usuaily expressed as 


C(s) G(s) ; 
Roa le Gls) 


Symmetrical two-dimensional, orthogonal systems 
containing no cross-coupling may generally be depicted 
by two such diagrams as, 


Ro(s) + Eo(s) C.(s) 

ee ee ee as = = 
R,(s) + E , (8) C , (s) 
a 


where 6 and ¢ denote the two orthogonal planes of ac- 
tion. Thus, the total transfer function may be described 
by two equations, 


(2) 


Cota serG 
eee o 
0 +G 
Ce iG 
pT es os 
$ Tats, 
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To depict linear cross-coupling in a more general 
fashion for a two-dimensional system, a general diagram 
for the total system may be represented as: 


where G(s) is now composed of G,(s) and G2(s) operat- 
ing in series. The signal present at the output of G,(s) 
is summed with a signal from the orthogonal channel. 
The summed signal is the input to G2(s), and is also 
operated on by /(s) to represent the signal cross-coupled 
between the two orthogonal systems. Algebraically, the 
following relations exist: 


Es = Ra — Co (5) 
Fe = E,Gi + FyJ (6) 
Co = Else (7) 
Hg ge Cy (8) 
Fy, = E,G, + FoJ (9) 
Ce = F Go. (10) 


These six equations may be solved to eliminate Fy, F4, 
Es, and E,. This solution provides the total transfer 
function for the system in two equations, 


2 : pit+e+ (11) 
Rate Gea J? R; 

G R 
eZ jite+ os]. (12) 
Rainey ad Ry 


These equations may be expressed in matrix form as, 


| E +G oi (fel G mies 
es at | eee fete tee i ted G)t sf? 
Should Rs = Ry, (11) and (12) reduce to 
ee eG 
ides tay (15) 
Kee taG 
Or, in matrix form, 
ba > Ba G (16) 
(ay RJ1+G-J 
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Should Rs= — Ry, (11) and (12) reduce to 


Co G 
= —— (17) 
Relea 
(On G 
= — . (18) 
Wi UNS Gs sea 
Or, in matrix form, 
Co Ro G 
= pene Comme | (19) 
Cy TEA 3) OKC dl 


The latter case is equivalent to the former case with J 
negative. The analysis of the systern performance thus 
becomes an analysis of the general case transfer func- 
tion which contains (1+G)?—J?=(1+G—J)(1+G+J) 
in the denominator, and the special case transfer func- 
tion which contains 1+G-—/J in the denominator, where 
J may be either positive or negative. 

Note that the location of the cross-coupling in the cir- 
cuits is immaterial. That is, the product, GiG2, always 
occurs and is equal to G, the forward transfer function, 
without cross-coupling. Inasmuch as G,; and G» do not 
appear independently, but only as products reducible 
to G, it will not be necessary to determine G,; and G» 
individually. 

Note further that in transfer functions of this nature, 
the reference term, Ry or Ry, is basically a normalizing 
term for the plane of action of the equation. Thus, any 
cross-coupling term generating from another input must 
occur in a normalized form (z.e., Ry effects in the Ro- 
plane appear as R,/R»). 

In spite of the generalized nature of the algebraic 
model thus far developed, there are some inherent limi- 
tations that bear mentioning. First, the model requires 
that each system component be describable in perform- 
ance by linear differential equations in real time. Second, 
the system has been assumed to be identical in both 
operating planes, and the cross-coupling transfer func- 
tion J has been assumed to be identical also. These as- 
sumptions have been made to simplify the treatment 
and the examples, and are considered justifiable on the 
basis that the majority of two-dimensional systems in 
existence are designed to be identical in both planes of 
action. Further, it is believed that the treatment pro- 
vided in this paper can provide guide lines and tech- 
niques to be followed in the analysis of non-symmetrical 
two-dimensional systems. 


III. STaBriLity 


In the analysis of cross-coupling effects on the per- 
formance of two-dimensional, orthogonal, feedback con- 
trol systems consideration must be given to the stability 
of the total system. This must be done in order to estab- 
lish whether the system is or is not stable, whether or 
not gain and phase margins exist, and to provide a basis 
for manipulating the response of the system so as to 
secure acceptable over-all response characteristics in 
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the presence of cross-coupling which cannot be elimi- 
nated. It is the purpose of this section to present an ex- 
position of how such an analysis may be made when the 
cross-coupling transfer function may be expressed by a 
linear, differential equation. 

Returning to the algebraic model developed in Sec- 
tion II, consider that it is a representation of an auto- 
matically tracking radar system in which the compo- 
nent and cross-coupling transfer functions may be fre- 
quency-sensitive. The general expression (13) may be 
investigated for system stability. 


eI pees uf i G 
Comme tee Git GS? 
To determine the stability of the general system rep- 
resented by (13), it is necessary first to establish that 
the zeros of the denominator are the zeros of (1+G)?— J? 


(z.e., no poles are contributed by the numerator). This 
may be done by letting 


(13) 


1 ng po. 
KG) d(J) 
Then 
n(G) n(G) n(G)n(J) 
Rp (rea eR eens 
on d(G) | cea d(G)d(J) (20) 


+56 - Gal 


This expression may be resolved to 


Ron(G)[d(G) +n(G) ||d/) }? + Ron(G)nJ)d(G)dV) 
Ce =— - (21) 
[4(G) +-n(G)]}?[4(J) ?— [nV P[a@) |? 


The zeros of the denominator of this expression are thus 
the zeros of the denominator of the general cross-coupled 
system transfer function (13), and the stability of the 
general system may be analyzed by studying the zeros 
of (1-+G)?—J?. 

It is necessary next, following Nyquist’s stability 
criterion, to determine that none of the zeros of 
(1+G)?—J? lie in the positive half of the complex s- 
plane (condition for stability). Although this may be 
interpreted as stating that an s=jw (Nyquist) plot of 
2G+G?—J? must follow the usual Nyquist rules 
regarding encirclement of the —1 point, the appli- 
cation of graphical techniques to determine a plot of 
2G+G*—J*? doesnot appear to be nearly so con- 
venient as it is to factor the expression (1+G)?—J? into 
(1+G—J)(1+G+J). For the factored form of the ex- 
pression to be equal to zero, one of the two factors must 
be equal to zero. These two possibilities are the same 
condition considered in Section II as special cases of the 
general system [(16) and (19) |. Therefore, stability con- 
ditions may be stipulated for these special cases, and 
these conditions should suffice to determine the stability 
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of the general system. That is, the system will be stable 
under all conditions if it is analyzed as a special case of 
the general system, and is determined to be stable for 
both positive and negative cross-coupling (+/J). 

In order to determine the stability of such a system 
it is necessary, following Nyquist’s stability criterion, 
to determine that none of the zeros of the denominator, 
1+G(s)+J(s), lie in the positive half of the complex 
s-plane (condition for stability), z.e., when 


Ie G(s)t-e 
(s — r1)(s — ro) + + + (8 — Tn) =0, (22) 


a iG ee 


it is necessary for stability that none of the roots, 
11, 72, 73) ° * * y Yn, have positive real parts. There is no 
such restriction on the other roots, 7a, 7%, %c, °° * 5 Ym: 
The stability criterion may be interpreted for this 
special case system as requiring that all —1 values of 
G(s) +J(s) be located in the left (or negative real) half 
of the s-plane. 

If a complex-plane curve for G(s)-+J(s) can be con- 
structed with s =7w for all values of wfrom —» to+, 
then the Nyquist stability criterion may be applied. A 
vector between the (—1+ 70) point to the (G(s) +J(s)) 
plot may be observed to determine the net number of 
counterclockwise revolutions of the vector as w varies 
from —« to +, and then determining whether this 
net number of counterclockwise revolutions is equal to 
the number of poles of the open loop transfer function 
in the right half plane (for stability) or not (instability). 
Note that this is an abbreviated version of the Nyquist 
stability criteria which includes the assumption that 
the net number of counterclockwise revolutions will be 
determined as stated in the procedure outlined on p. 145 
of [7]. A zero of 1+G(s)+J(s) with positive real parts 
will produce a clockwise revolution of the vector, and a 
pole of 1+G(s)+/J(s) with positive real parts will pro- 
duce a counterclockwise revolution of the vector. It is 
usually assumed that G(s) and J(s) are by themselves 
stable, 2.€., 7a, 7, * * * » Ym, the poles of 1+G(s)+J(s) all 
lie in the negative half plane. For systems in which 
1+G(s)+J(s) contains poles in the positive half plane, 
it is necessary for stability that the number of counter- 
clockwise revolutions of the Nyquist vector equal the 
number of poles in the positive half plane. Note that 
the following two limitations are implicit in using the 
Nyquist stability criteria; 1) the system must be capable 
of being represented by a system of linear differential 
equations with constant coefficients, and 2) the limit of 
G(s), and of J(s), must approach a constant or zero as s 
approaches ». The latter limitation is inherently pres- 
ent physically in linear systems. Mathematically, it is 
required that for any physical system the transfer func- 
tion contain a higher power polynomial in its denomina- 
tor than in its numerator. 

In order to make a complex plane plot of the function 
G(s) +J(s) with s=jw, and allowing w to go from — « 
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to + «, the following procedure appears to be the least 
tedious available with presently known mathematical 
techniques: 

1) Make an s-plane plot of the poles and zeros for 
G(s). By graphical construction, take values of s =jw as 
w goes from — ~ to + ~, determine the magnitude and 
phase angle for G(jw), and plot these values on a com- 
plex plane. This will produce the typical Nyquist dia- 
gram, which is ordinarily analyzed for stability by de- 
termining enclosures of the —1 point, or the number of 
rotations of the vector between (—1+ 0) and the 
G(jw) curve as w goes from — © to +o. 

2) Repeat the foregoing procedure to obtain a simi- 
lar plot for J(jw) on the same complex plane. 

3) Indicate identical values of w at various points on 
both the G(jw) and J(jw) plots, thus producing a plot 
such as Fig. 1. 

4) Sum the two curves vectorially for several values 
of w, and make a new complex plane plot which will be 
a plot of G(jw) +J(jw) as w goes from — ~ to +. The 
plot may have the form of Fig. 2. 

This plot may then be analyzed relative to the 
—1-+0 point following the usual Nyquist method. 

The example used depicted a J(s) of the nature 


Jx 


a § = ) 23 
Bee eae tt oe 
and a G(s) of the nature 
Gr 
G(s) (24) 


Eee (hie 1). 


This combination does not appear to give rise to any 
serious stability problems. However, had J(s) been of 
the nature 


Jk 
= DS 
ve) $(taS + 1) (ros + 1)e7* oe 
or 
Tx(ras + 1)(res + 1) (tes + 1) 

— d 26 
ete, Ga 1) i?) 

with 


Ta, Tb, Te ee Td, Tey Tf; 


a serious stability problem might occur. 

Consider the individual plots of step 3) in the pro- 
cedure shown in Fig. 3 (next page). 

A summation of these two plots may yield a plot such 
as Fig. 4. Fig. 4 indicates how the mere addition of the 
cross-coupling transfer function may easily lead to an 
absolute instability which may or may not be elimi- 
nated by a simple gain decrease in G(s), the stability be- 
ing dependent also on the gain and phase characteris- 
tics of J(s). A similar summation should be made for the 
other special case obtained by taking J as negative in 
sign. Such a plot is depicted in Fig. 5. 
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a 
~I 


6S) = 
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T Jk 
J(S ) = ———————— 
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Fig. 1—Complex-plane plots of Gw) and J(jw). 
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Fig. 2—Complex-plane plot of G(jw) +J(jw); system not 
made unstable by J addition. 
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GK 
UY eee eee 
S$(T,S +1)(T2S41) 
J 
J(s) = ‘ 


S(TQS+1)(TpSt!) 


Fig. 3—Complex-plane plots of G(jw) and J(jw). 


It should be remembered that both special cases 
(+J) must be stable for the system to be stable, and 
that the use of cross-coupling to make one case more 
stable will generally tend to make the companion case 
less stable. The stability margin of the total system will 
usually not be greater than either special case (+/J) 
stability margin. 

To illustrate the analytical methods of this section 
further, consider a positional system, 


0.5(0.1s +1) 


(0. 5s +1) 


0.5(0.1s +1) 


(0. 5s + 1) 


BOs Hs 2.0 


G(s) = . 
(s+ 1)(2s+ 1) 35+ 1) 
es 10.0(0.2s + 1) (21) 
Ce IE PN ieee 8 
oe 0.5(0.1s + 1) (28) 


(0.55 + 1) 


Gy 


S(T, st!)(Tas+1) 


vk 
S(T, st1)(Tps +1) 


J(s)= 


Fig. 4—Complex-plane plot of G(jw) +J(jw): system 
made unstable by J addition. 
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vk 
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Fig. 5—Complex-plane plot of G(jw) —J(jw); system not 
made unstable by J subtraction. 
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10.0(0.25+1) 
(S#1)(2S+1)(3S+1) 


G(s) = 


Ja) = OS (OSH) j 
(0.5 St!) 


Fig. 6—Complex-plane plot of G(jw) and J(jw) for 
time constant example. 


G(jw) and J(jw) may be plotted as illustrated in Fig. 
6. Then, from these plots we may plot G(jw) +J(jw), 
Fig. 7, and G(jw) —J(jw), Fig. 8. These plots illustrate 
again, using actual time constant values, the manner 
in which cross-coupling may cause a stable system to 
become unstable. 

A demonstration of how cross-coupling may be of a 
beneficial nature in a system, and therefore purposely 
introduced, may be obtained by considering a system 
whose open loop, one axis transfer function may be 


Pa Gr(rT35 i 1) (745 = 1) 
s(715 —- 1) (ros oo 1) (755 + 1) 


G(s) (29) 


with 
‘il. ae it) ei ee itp 
and such that a complex plane plot will appear as de- 
picted in Fig. 9 (next page). Such a system may be im- 
mediately judged as unstable. However, consider intro- 
troducing cross-coupling of the nature, 
Jx(ras + 1)(re5s + 1) (7s + 1) 


s 30 
eee yg tt) (30) 


with 
Ue Ty) ee Uy Lee 
and Jx>1, and such that a complex plane plot will ap- 


pear as depicted in Fig. 9. Making the graphical sum- 
mation of G(jw) + J(jw) should then yield the complex 
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aleve 10.0(0.2S +1) 
(S+1)(2S+1)(3S 41) -j 
0.5 (0.1S +1) 
J(s) = —————_——_ 
(0.5S+1) 


Fig. 7—Complex-plane plot of G(jw#) +J (jw) for 
time constant example. 
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Cieie 10.0(0.2S+1) 
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w(s) =——————— 
(0.5S+1) =] 


Fig. 8—Complex-plane plot of G(jw) —J(jw) for 
time constant example. 
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Fig. 9—Complex-plane plot of unstable G(jw) and J(jw) 
introduced to stabilize system. 


plane plots depicted in Figs. 10 and 11; and by com- 
parison with Fig. 9 it may be observed that the unstable 
system has been made stable by the introduction of 
cross-coupling. Thus, cross-coupling may be usefully 
introduced for the purpose of making an unstable sys- 
tem stable (although consideration must be given to 
the seriousness of adverse effects such as error vector 
distortion). 


IV. CONCLUSIONS 


The stability of cross-coupled systems may be inves- 
tigated by graphical techniques following the procedure 
outlined in Section III. The total system stability is de- 
pendent on the forward transfer function and the cross- 
coupling transfer function. The presence of cross- 
coupling in an otherwise stable two-dimensional system 
may cause the system to become unstable. Unstable 
two-dimensional systems may be made stable by the 
introduction of suitable cross-coupling transfer func- 
tions to the system, as illustrated in Section III. 

The work of other authors, particularly [1]|-[6], has 
been added to by this paper in that a general treatment 
is now provided for the analysis of the stability of any 
cross-coupled system in which the cross-coupling trans- 
fer function may be expressed by a linear differential 
equation, and the analysis may be accomplished by 
complex plane methods. 
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Stabilization of Linear Multivariable Feedback 
Control Systems” 


E. V. BOHN, AssoclaTE MEMBER, IRE 


Summary—Stabilization techniques for single variable feedback 
control systems are well known. For multivariable feedback control 
systems, comparable methods are still in the initial stage of develop- 
ment. A method of system stabilization is discussed, which reduces 
the problem to essentially that of a single variable system. This 
method uses a compensating matrix which is either the inverse sys- 
tem- or transposed system-matrix. An example of this method of 
system stabilization is given and its usefulness discussed. 


INTRODUCTION 


TABILIZATION techniques for single variable 
S feedback control—or servo—systems have been 

brought to a high degree of development. For mul- 
tivariable servo systems, the development of compar- 
able techniques is still in its initial stages. It is of interest 
to give a brief discussion of the relevant material avail- 
able in the literature. Consider the linear multivariable 
servo system shown in Fig. 1. X, Y and V are single 
column matrices of » elements each where X is the 
input, Y is the controlled output, and V is the external 
disturbance. A, C, G and H are square matrices of n? 
elements each. A represents the system; C is the com- 
pensation matrix; G represents the high gain power 
amplifier stages; and H is a feedback matrix. The ele- 
ments of these matrices are the Laplace transforms of 
the respective time functions. These matrices are re- 
lated by the following equations: 


U =C(X — 2), (1) 
W=UG+Y, (2) 
Y = AW, (3) 
Z = BY. (4) 


By a process of elimination it is easy to obtain 
(I + AGCH)Y = (AGC)X + AV, (5) 


where J is the unit matrix. 

A system of this type, where G and H are diagonal, 
has been extensively discussed in the literature as an 
example of noninteracting controls.'~* The condition for 
noninteraction is X =Z. From the servo system point of 
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Inc., New York, N. Y., pp. 53-69; 1954. re: ars, 

2R, Kavanagh, “Noninteracting controls in linear multivaria e 
systems,” Trans. AIEE, vol. 76, pt. 2, pp. 95-100; May, 1957. 
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Fig. 1—Multivariable servo system. 


view, this is achieved by having a high loop gain. The 
system is then relatively insensitive to both the external 
disturbance matrix V and to the interaction of signals. 
In this sense, every multivariable servo system is non- 
interacting. However, this is not the approach taken in 
the references cited. To understand this situation, con- 
sider the problem of determining a compensating matrix 
C in order to obtain a stable system with a high gain in 
the elements of G. One possibility is to choose 


Cit DAR® (6) 
where D is diagonal. For this choice of C we have 


AGCH = DGH, (7) 


which is diagonal. Eq. (5) can now be written as 
DikGur (AV) 
ts Ah fae hee oe 
1+ DixGer Hen 1+ DicGurH ex ’ 
| Sige Wl oan etait 5) 


Vx 


This can be simplified to 


VY, = RipXn — (RerHiz — 1)( AV); (9) 
where 
DyrGir 


Ki Se ae 
1 Du-Guls 


(10) 


This particular choice of C diagonalizes the open loop 
transfer matrix AGCH. Eq. (9) can now be treated 
as nm independent servo systems, and conventional 
single variable techniques can be used to determine 
suitable compensating elements for Dx, or Gzx. These 
single variable techniques are the extensively developed 
methods based on the root locus, Nyquist plots, and 
phase and gain margins. 

The result (9) is given in the cited references as the 
condition for noninteraction. This terminology is con- 
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fusing since every high gain multivariable servo system 
is noninteracting. H. Freeman, in the discussion of 
Kavanagh's paper,” has criticized the interpretation of 
noninteraction from a different point of view. 

In a separate paper, Freeman‘ discusses the synthesis 
problem and introduces the concept of “independent 
output restoration.” This terminology is unnecessary 
since the open loop transfer matrix is diagonalized ex- 
actly as is done in the system just discussed. The differ- 
ence is simply that, in Freeman’s case, His not diagonal. 
The synthesis and related problems are also discussed 
by Freeman® and Kavanagh.® Both these latter papers 
rely on having determined, by some means, the desired 
closed loop transfer matrix K. A compensation matrix 
C is then determined to realize this K. This approach 
has been criticized by F. Nesline in the discussion of 
Kavanagh’s paper.® Before a suitable choice of K can be 
made, the dynamics of the closed loop system must be 
identified and approximated. It is necessary to investi- 
gate the response of the fixed elements of the system 
with respect to some criterion. A simple criterion is that 
of having a high loop gain consistent with a stable dy- 
namic response. The question of stability and response 
of multivariable servo systems has received little atten- 
tion. Very special cases are discussed by Golomb and 
Usdin’ and Pelegrin.® 


MULTIVARIABLE SERVO SYSTEMS 


To avoid undue complexity in discussing the stability 
problem, we will set V=0 and H=/. The system matrix 


1S 
Ain i ie Ain 


(11) 


ue ‘ wtAgs 


where the A;;’s are considered to be fixed transfer func- 
tions. G will be taken to be diagonal. The simplifying 
assumption is made that all power amplifier stages are 
identical and have a transfer function G(s). This in- 
volves no restrictions since, through suitable compen- 
sation networks, this can always be achieved. Hence 


CSCC RE (12) 


The uncompensated multivariable servo system is 
shown in Fig. 2. The matrix equations for this system 
are 


4H. Freeman, “A synthesis method of multipole control systems,” 
Trans. AIEE, vol. 76, pt. 2, pp. 28-31; March, 1957. 

5 H. Freeman, “Stability and physical realizability considerations 
in the synthesis of multipole control systems,” Trans. AIEE, vol. 
77, pt. 2, pp. 1-5; March, 1958. 

6 R. Kavanagh, “Multivariable control system synthesis,” Trans. 
AIEE, vol. 77, pt. 2, pp. 425-429; November, 1958. 

7M. Golomb and E. Usdin, “A theory of multidimensional servo 
systems,” J. Franklin Inst., vol. 253, pp. 28-57; January, 1952. 

8 J. Gill, M. Pelegrin, and P. Decaulne, “Feedback Control Sys- 
tems,” McGraw-Hill Book Co., Inc., New York, N. Y., pp. 370-373; 
1959. 
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V = TARY. (13) 
H=V-—-Y=VJ — AX, (14) 
E 
= G()b = =, (15) 
d 
where we have set 
r (16) 
Pare, 
From (13)—(16), it is easy to obtain 
(A AD na Va (17) 


Eq. (17) will be used in the discussion of system sta- 
bility. 
SYSTEM STABILITY 


The stability of the system can be determined by 
means of the Nyquist criterion. Let 


Ae(S) 5 R=1j:.-<,” 


be the eigenvalues of the A-matrix, 7.e., the values of \ 
which satisfy the equation 


A— d-I| =0. (18) 
The system is stable if there are no values of s in the 
right half s-plane which are roots of (18), since these 
are the poles of the closed loop system. To determine 
whether this is the case, consider the function R,(s) 
defined by 


1 


Ri(s) = uo) 


G(s): (19) 


From (16), (18), and (19) it is easily seen that the roots 
of (18) are the zeros of 


Ri(s) = 0. (20) 
The condition for stability is that the zeros of (20) must 
all lie in the left half s-plane. This can be determined 
by considering the Nyquist plot of R;(jw) and applica- 
tion of 


N=Z-P, (21) 
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where P is the number of poles and Z the number of 
zeros of R;(s) in the right half s-plane respectively. N is 
the number of positive revolutions of the radius vector 
R,(jw). If Z=0, the system is stable; if Z>0, there are 
values of s in the right half s-plane which are roots of 
(18), and the system is unstable. 


SYSTEM STABILIZATION 


A general theory of the stabilization of multivariable 
servo systems would be very complicated, and its use- 
fulness for design purposes limited, due to the compu- 
tational difficulties involved. A method of multivari- 
able servo system stabilization is discussed which re- 
duces the problem to essentially that of a single variable 
system. It is shown how this objective can be accom- 
plished byasuitable choiceof the compensating matrix C. 

The system already discussed with reference to Fig. 1 
and (8) indicates one possibility. Eq. (6) shows that the 
compensating matrix is essentially the inverse of the 
system matrix. This coice of C leads to » independent 
systems. The stability and response of each of these can 
be investigated by conventional single variable tech- 
niques and suitable compensating elements determined. 
This method appears to be very general. Some of its pos- 
sible disadvantages in a practical system will be dis- 
cussed later. The choice C=A7 for the compensation 
matrix will be called type I compensation. 

Consider now the problem of stabilizing the system 
discussed with reference to Fig. 2 and (16) and (17). 
The system eigenvalues \;(s) may lie anywhere in the 
s-plane. To ensure system stability and reduce the com- 
plexity of the stabilization problem, the C-matrix will 
be chosen so as to localize the system eigenvalues to 
specified regions of the s-plane. To see how this can be 
accomplished consider the augmented system 


(B-r-(")= (2), (22) 
W V 
where 
eee). @ 
This is equivalent to the two matrix equations, 
CW = (1—A)X, (24) 
AX =)W-Y. (25) 


There are two possible interpretations of this augmented 
system. One is similar to Fig. 2, with V, X and A being 


replaced by 
(a) ( - 
Vv)’ W 


and B, respectively. This representation is based on (22). 
Eqs. (24) and (25) lead to the representation in Fig. 3. 
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Fig. 3—Augmented multivariable servo system. 


[It is important to note that d is defined by (16).] The 
use of an augmented matrix of this type was introduced 
by Goldberg and Brown! to stabilize the solution of 
linear equations on an analog computer. In this appli- 
cation, the system matrix A contained only real num- 
bers, and C was chosen to equal A’, the transpose of A. 

Eqs. (24) and (25) can be simplified by eliminating 
W. This leads to 


(CA — y1)X = CD, (26) 

where 
ray V Ghee oO): (27) 
Consider now the eigenvalues y;(kR=1,---,) of 


the matrix CA. If ; is an eigenvalue of the augmented 
matrix B, (27) gives 


Ye = Ae(L — Ax). (28) 
This can be solved for A; to give 
Ae = 2(L & iV 47x — 1). (29) 


To see how C can restrict the location of the eigen- 
values of B in the s-plane, suppose that CA is a positive 
hermitian matrix. Then y; is real and positive (see Ap- 
pendix III). The eigenvalue locus for all augmented 
systems is then located on the two straight lines. 


VAN 
BIR BIR 


H piers O<” a 


(Ad + jar — 1) ies 


These two lines will be called the real locus. 

Since the elements of the C-matrix must consist of 
physically realizable transfer functions, the question 
immediately arises as to how this is to be achieved in 
addition to ensuring that CA be hermitian. A simple 
possible choice for C is 


C= hAM, (31) 


where his a real constant >0. This is type I compensa- 
tion. In this case, CA =hI. Hence yi =h and 


he = 3(1 + V1 — 4h). (32) 


9 E. Goldberg and G. Brown, “An electronic simultaneous equa- 
tion solver,” J, Appl. Phys., vol. 19, pp. 338-345; April, 1948. 
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Let \1 and A, be the two values of the right hand side of 
(32). The two points, —1/A, and —1/), in the s-plane 
will be called critical points corresponding to the point 
—1 in the one variabie case. The conventional single 
variable Nyquist criterion is now applicable. This 
means that, provided P =0 [see (21) |, the encirclement 
of the critical points by the locus G(jw) is to be avoided 
by a suitable gain and phase margin. 

A second possible choice is C=A*, where A* is the 
adjoint of A (see Appendix III). However, the elements 
of A* are not generally physically realizable. The con- 
ditions can be partially met by simply choosing C=<A?. 
For s real A* = A7™ and the A; (s) locus will lie on the real 
locus. For s=jw, the Nyquist locus \;(jw) will be sym- 
metrical about the real locus as shown in Fig. 4 (for 
simplicity only half of the locus is drawn). Hence, this 
choice of C also restricts the eigenvalue locus. This type 
of compensation will be called type II. 

The stability can once more be determined by appli- 
cation of (21). For example, suppose P =0 and the Ny- 
quist plots of G(jw) and —1/)\;(jw) are as shown in 
Fig. 4. A locus must be plotted for each R=1,---, 7. 
Z=N is determined in each case from the number of 
positive revolutions of the radius vector R,(jw). A 
simple approximation can be obtained by considering 
the largest value max of all the \;’s on the real locus. 
—1/Nmax can then be taken as the critical point and 
single variable techniques may be usefully applied. 


PHYSICAL REALIZABILITY 


With type I compensation, C is physically realizable 
if the elements of A~! are physically realizable as trans- 
fer functions. This is the case if A! has no poles (and 
hence if A has no zeros) in the right half s-plane. The 
C-matrix with type II compensation is always physically 
realizable. 


CONCLUSION 


By suitable choice of the transfer function of the 
power amplifier GOS), type I compensation will evi- 
dently always lead to a stable system, provided A has 
no zeros in the right half s-plane. The dynamic response 
may be somewhat sluggish due to the multiplicity of 
equal roots. This can be seen from (32). The augmented 
system has eigenvalues each equal to \y, and n eigen- 
values each equal to de, respectively. The zeros of (20), 
which are the poles of the closed loop system, will have 
a corresponding multiplicity. A practical disadvantage 
of this method is that it evidently depends on having 
the actual system closely approximated by a linear one 
over the entire range of operation. 

Type II compensation may not always lead to a 
stable system. However, it has the practical advantage 
in that the C-matrix may consist of physical analogs or 
models of the system. The method can then be applied 
to nonlinear multivariable servo systems which can be 
linearized by describing function methods. 
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Fig. 4—)-plane. 


APPENDIX I| 


As an example, consider a system whose matrix is 


A= (s 3) (33) 
1 ws 
where 
OS eae . (34) 
1] 
For the power amplifier, assume 
coe es (35) 
s(1 + sr) 


The performance of this multivariable system will be 
judged by comparing it with a reference system. This 
system is chosen to be a conventional single variable 
servo system having the same power amplifier. For the 
reference system, the open loop transfer function is then 
given by (35). If a relative damping ratio of ¢=0.6 is 
chosen for the poles of the closed loop reference system, 
it is found that 


Kyr = 0.945, (36) 
The velocity error for a unit ramp input is 
e1(©) = 1.067. (37) 


To determine the stability of the unmodified multi- 
variable servo system, the system eigenvalues, \; and Xa», 
must first be found from (18). These are 


M=V4+; 


hase VEG (38) 
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The poles of the closed loop system are the roots of 


1 s(1 + sr) 1 
MI = us ae. 
1+ sT 


Tie (39) 


G(s) I 
and a similar equation with A». By substituting sT=w, 
t/T=a, and K,T=b, the following cubic equation is 
obtained: 


u(i+ auj(i+ u)+7b4+u)+6=0. (40) 


Assume now that a<1. One root of (40) is then approxi- 
mately —1/a, and the other two roots are solutions of 
the quadratic 


w+ it+jbu+ i+ 7)b = 0. (41) 


The gain K, is to be chosen so that the relative damping 
ratio of the roots of the system will be 20.6. Solving 
(41) under this constraint leads to the following results: 


b = 0.29 

a=) —.0.28 — 90.47; 

te = — 0.72 + 70.18, 
1 

4u,=—-—-: (42) 
a 


The eigenvalue dy» will evidently result in the conjugate 
complex of these roots. 
Since K,7T =), we have 


K,r = 0.29-a. (43) 


The gain is extremely small compared to the reference 
system, and the performance of the system will be poor. 
This can also be seen from the Nyquist plot of the un- 
compensated system in Fig. 5. It is evident that a con- 
siderable reduction in gain is required to ensure that 
N=0. 


Type I Compensation 

The system is augmented with the compensation 
matrix C=1A—. The augmented matrix B has four 
eigenvalues A;.(k=1, 2, 3, 4) [see (22) ]. From (31) and 
(32), we have h=} and \;,=}. The poles of the closed 
loop system are the roots of 


+r. = 0. (44) 
G(s) 
Substituting G(s) and using [=0.6 gives 
Kar = 1.89. (45) 


The velocity error matrix can be obtained from (68), in 


Appendix II, where (AC) +=) *=4I. Hence, 


(46) 


(0) = —U = 2.1270, 


K, 


where U is a single column unit matrix. 


s — 
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Fig. 5—)-plane, uncompensated system. 


The augmented system has the same relative damping 
ratio and double the velocity error of the reference sys- 
tem. System performance can be further improved by 
shaping the G-locus. 

It is of interest to determine the elements of the C- 
matrix for type I compensation. 

In this case, 


Ain Aa 
A A 
AG’ = A= ‘ (47) 
Ai Arg 
A A 


where the A,;,’s are the cofactors of A and A= | A | . We 
have 


Auta) VY _ Ave 

Aton tay oA 

A 4 A 

ike epee aeg (48) 
A ae 


where Y=1/(1+u) and w=sT. 
Since the elements of C are to be realized by means of 
transfer functions, one possibility is to set [see Fig. 6(a) | 


i} i 
Au as nee Goi Mae hs Se tiebeg - (49) 
A Zi + Zo Ze U? + U + 1 
1+— 

1 U + 1 

Arbitrarily choosing Z,=1 leads to [see Fig. 6(b) | 
1 
Li i (50) 
aes 1+u 
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ot | 
7 I 
(b) (c) 


Fig. 6—Elements of the compensation matrix X. 


Similarly, 


Aoi 1 1 
‘ (S1) 


LZ. 
1+ — 1+ 
Z\ (1 + a)? 


Arbitrarily choosing Z,;=1/(1+u) leads to [see Fig. 
6(c) | 


A=1+4. (52) 


Type II Compensation 
In this case, C= AT and hence 
1+ Yy 0 


CA = Ard ~ ( 
0 As we a 


) Se( ds V2) 0 (53) 
The eigenvalues for CA are yi=Yy2=14+ Y*. The eigen- 


values \;(s) for the augmented matrix B are determined 
from (29). Since 


Ky b 

s(1 + sr) . u(1 + au) 
and since a has been assumed small, wu will be large in 
the region of interest (w1/a). This means that the 
eigenvalue locus A;(jw) need not be completely deter- 
mined, and a good approximation is obtained by 
simply considering the value \;(j ©). Since Y(j~) =0, 
¥1=7Y2=1, and (29) gives 


(jo) = $1 + jV3], 


G(s) = (54) 


(55) 


the two points 


1 
= ey LEENA (56) 


can then be considered as critical points, and conven- 
tional single variable techniques applied. 


APPENDIX II 


Error coefficients can be used as a performance index 
for a stable system. For a one variable system the error 
is 


SV SG he Gal Vix: oY 
a Seah ete Gala Ve (57) 
The velocity error for a unit ramp input is 
ej() = lim ski(s). (58) 
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For zero position error we must have 


1 


lim = 0; (59) 
s—0 G(s) 
hence 
eG 3) 1 
€:(%©) = lim inf oS , (60) 
s—0 a Ke 


where K, is the velocity error coefficient. 

For the augmented multivariable servo system, the 
error matrix is readily determined from Fig. 3, and is 
given by 


B= — dW. (61) 


Eqs. (24) and (27) can be used to eliminate W and X. 
The result is then 


E> — CX, (02) 
From (26) we obtain 
KEa(CA = ye eCie (63) 
Using the identity 
(CA =~ yDOUOC = = yA GS) ASCH 
=! (i ay Ca etd (64) 


and eliminating X from (62) and (63), it is easy to 
obtain 


E=—yCtANT — yOrd ty, (65) 
For zero position error we must have 
lim a — yC-1A71 = 0: (66) 
Hence, for unit ramp inputs, 
e(co) = lim on — ae (67) 


where €(~) is the velocity error matrix and U is a unit © 
column matrix. Eqs. (65) and (67) correspond to (57) 
and (60) in the one variable case. Eq. (67) may be fur- 
ther simplified. From (16) 


r 1 1 
—>— as s— 0. 
5 SG(S) 5 Ky 
Hence A-0 and y—A. Then 
y d 1 
Ss Ss IK 


The velocity error matrix is then 


MOE — lim inf (AC)U. (68) 


v s>0 
By a similar analysis it may be shown that (68) is also 
valid for the uncompensated system, provided C is re- 
placed by J. 
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APPENDIX III 


A few of the results from linear matrix—or vector— 
algebra are included for reference purposes: 


The scalar product of two vectors u and v can be de- 
fined by 


(u, v) = Uy*01 + pn =e Un* Un = (v, u)*, (69) 


where 4, - ++, u, are the elements of uw, and u;* is the 
conjugate of uz. Then 


(u,u) = | m4 


24 ..-4 | u, 


2> 0. (70) 


The adjoint B* of a matrix B is defined by 
(u, Bo) = (B*u, 2). (71) 


This definition is that used by Laning and Battin.’ 
It is easily verified that B*(s) = B?(s*), where B® is the 
transpose of B. If B*=B, the matrix is hermitian: 


10 J. Laning and R. Battin, “Random Processes in Automatic 
Control,” McGraw-Hill Book Co., Inc., New York, N. Y., pp. 332- 
333; 1956. 


Bohn: Stabilization of Linear Multivariable Feedback Control Systems 


Consider the eigenvalue equation 
Bu = Xx; 
scalar mutliplying by x gives 
(Gauls) NCA ale 
If B is hermitian, (x, Bx) is real, since 
(CODD a= OD ee ese Ys eye 07 eee 


Hence 


$B, JES 

(x, Bx) So 

(x; a) 

We can restrict ourselves to values 20, since, if A 


we can replace B by —B. 
If 


pS Med. 
then 
(Bury = (A* Ax a) = (Ag An) =i xB) 210. 


Hence B is hermitian and ) is real and >0. 


a2, 


(74) 


S0, 
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Comments on “Optimization Based 
on a Square Error Criterion”* 


The main results in the paper by Murphy 
and Bold! appear to be incorrect. The first 
one is (46) of the paper, which gives the 
optimum system frequency function as 


. \ _ Vucgr( Jor, jo’) 
(Gs se eee ed Fe 1 
pt(jo) WVor( jo, jw’) ( ) 


This is the case where the question of phys- 
ical realizability is disregarded. The func- 
tions Vycgr and V,, are the double Fourier 
transforms of the weighted correlation func- 
tions corresponding to the desired output 
ca(t) and the input r(é). Since their result is 
based on minimizing the weighted square 
of the error 
1 


: F 
B= lim = J Woe (dt, (2) 


then, for the case where W(t) =1, their re- 
sult should be reducible to the classical 
Wiener result for pure smoothing, 


Steg) Sow) 
Sr(jo) Sew) + Sw(@) 


Gopt (jw) =! (3 ) 


Consider the case of pure smoothing of 
a signal in the presence of white noise. Let 
the autocorrelation function of the signal 
input be 
g(t) =e"! =el*I, forallr,*,andy. (4) 


The corresponding signal power spectrum is 


0 c-) 
Ss(w) = a ee id + df Cen dt 
-0 0 


2 


eae ©) 


The double Fourier transform of e7l2-v! 
with respect to x and y, however, does not 
exist. We have 


eo 
f elt wleI6% J, 


y eo 
= er Ve I dye + i ete 1rd y 
y 


(6) 


But 
race is) es i ei g—vHq (7) 
rr\J, J Bei at é iy 


does not exist. Hence, (46) is meaningless for 
this example. 

The root of the trouble seems to be in the 
step taken by the authors in going from (39) 
to (40). Their (39) is 


* Received by the PGAC, April 11, 1960. 

1G, J. Murphy and N. T. Bold, “Optimization 
based On‘a square-error criterion with an arbitrary 
weighting function,” IRE TRANs. ON AUTOMATIC 
ContROL, vol. AC-5, pp, 24-30; January, 1960, 


i} 86) [Cee a 


mF bw, (T re aie y) dy = Puegr(T, ae x) 


= 0, (8) 


T=y=0 


sea dwegr (Ys Tae x) 


for « greater than zero. Note that the left- 
hand side of the above equation is inde- 
pendent of 7 and y. The authors claim that 
taking the double Fourier transform of both 
sides of their (39) yields their (40), 


J ye ‘| an f s (9) [ur(y—2, >=) 


al en Ta y) |dyd-ydr 
_ f call CF buce(T, y — x) 


ey 1) |ayart (0) 


T=y=0 


This is, of course, incorrect. Since the left- 
hand side of (8) is independent of 7 and y, 
the Fourier transform, with respect to 7 and 
y or both does not even exist. The rest of 
their equations leading to (46) are, there- 
fore, questionable. 

The principal contribution of the paper 
seems to be the derivation of their (39), 
which may be simplified further as 


J 80 )lée(—2, ~ lay = dene, —2), 


forx>0. (10) 


For physical realizability, g(y) must be zero 
for y negative, so that (10) becomes 


if 2(9) [bere(—a, —9) dy = ducye(0, —2), 


forse) Ose Gili) 


This is, of course, a generalized Wiener- 
Hopf integral equation. The weighted cor- 
relation functions are 


urr(ti, t2) 
1 T 
= jim J Wore + aret ai (12) 


and 
Pwegr G ? ta) 


1 iy 
= jim > JW Ocalt + trie + tae (13) 


Note that for 
Wii for 
(11) reduces to 


SES Fe Seco! | (GIES) 


if “Lodbete — y)dy = deyr(—x), (15) 


which is the usual Wiener-Hopf integral 
equation. 

Discussions with Dr. M. E. Van Valken- 
burg are greatly appreciated. 

jmBeCRruzeiR, 

Dept. of Elec. Engrg. 

and Coordinated Sci. Lab. 

University of Illinois 

Urbana, IIl. 
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Author's Comment? 


I am pleased to note the interest of Cruz 
and his colleagues in the paper under dis- 
cussion, and I am sorry to learn of their 
difficulty in understanding the material in 
the paper. However, I wish to point out that 
this difficulty was caused not by an error in 
the work described in the paper,! but rather 
by the use of a somewhat uncommon nota- 
tion, which was employed, in spite of its 
shortcoming, for the purpose of reducing the 
complexity of (40) in the paper. 

An alternative to the approach followed 
in the paper is to rewrite (39) in the form 


[ 52°:} f e0)lbmnty — 79) 


+ durr(t — x, — y) lay a Pwegr(T 7 — 2) 


= deel — 2) | = 0, «>0, 


T=y=0 


(39a) 


provided that the double Fourier transform 
denoted by S2{ } exists. An equivalent 
form is 


Salas 
des (hog 
47? —«0 Y —o —~ —« 


vei f 8(9) [durex — @, 7 — y) 
+ durn(r — x, 7 — 9) |dydydr 
= f eet { C791 bye an(7, y-x) 


=F Pweg (Ys T — x) |dydr dwdw' = 0, 


x>0, (39b) 
provided again that the double Fourier 
transform exists. Following the steps out- 
lined on page 28 of the paper leads from 
(39b) to (43) of the paper. Then, as is shown 
in the paper, the optimum frequency func- 
tion is found to be given by (46) of the 
paper. 

The example presented in (4)-(7) of 
Mr. Cruz’s discussion is based on the as- 
sumption that W(¢)=1. In this case, the 
two-space correlation functions reduce to 
the ordinary correlation functions, the prob- 
lem reduces to the usual problem, and there 
is no need to consider double Fourier trans- 
forms at all. In any event, it is clearly 
stated in the paper that (40) follows from 
(39) only if the integral is absolutely con- 
vergent. Obviously, Mr. Cruz chose an 
example in which this condition is not met. 
Hence his criticism is not justified. 

It should be noted that merely by choos- 
ing a sinusoidal form for r(é) one could 
similarly “show” that the well-known solu- 
tion to the ordinary Wiener-Hopf integral 
equation obtained through the usé of either 
the Bode-Shannon method or the Wiener- 


2 Received by the PGAC, April 25, 1960. 
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Kolmogoroff method is “meaningless.” 
Clearly, in this example as well as in the 
argument presented by Mr. Cruz, the fault 
lies in the choice of an input function for 
which the method is not rigorously appli- 
cable. 

: Alternatively, as is sometimes done in 
using the Wiener-Kolmogoroff method, one 
can resort to the concept of impulse func- 
tions, writing 


1 -) 
es i iid = uy(r). (66) 
Baie 


Then (7) of Mr. Cruz’s example becomes 


Wiss(iw, jor’) = i io ff cwee dy (67) 
— 4n , 8 
Fay aes (68) 
For the white noise one has 
gnn(T) = Uo(r) (69) 
dnn(x — y) = t.o(x% — y) (70) 
and 
Wrn(jo) = af ety (x — y)dx (71) 
= (ee (72) 
Then 
TWrn(Jo) = Vunn( jo, jw’) (73) 
= f e-7° Ve iV dy (74) 
= 27(w +’). (75) 
It follows from (49) that 
Gopt (jo) 
Galje) [7 we +2" | 
AGT) ite Ug(w + w 
= (76) 
a nletualy Bonita) 
1 xa wo? Ug\@ @ 0 
2, 
1 2 
= a Guljw). (77) 
1 
1+ wo? a 


Using the usual (Wiener) approach, one 
has 


P Ga (jo) Ves G w) 
je) = 78 
Se G ) Ws (jw) == Win(Jw) ( ) 
2 
1+ w? 
= 1 — Gulf). (19 
1+ o? z¥ 


Obviously, the result obtained in (77) is 
identical to that obtained in (79). Thus, it 
has been shown that, contrary to Mr. Cruz’s 
statement for the example chosen by him, 
the result obtained by application of the 
method presented in the paper is “reducible 
to the classical Wiener result.” 

I fail to see any point in the inclusion of 
(12) and (13) of the discussion, since these 
equations are identical to (26) and (25); 
respectively, of the paper. 
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The statement that “for 


W@) =1, for —~ <i< 0 


(11) reduces to 
f &(V) ber(x — y)dy = deyr(—2x), 
0 


which is the usual Wiener-Hopf integral 

equation,” was implied by (27) through 
(30) and the associated text in the paper. 

G. J. MurpHy 

Dept. of Elec. Engrg. 

Northwestern University 

Evanston, III. 


Nonlinear Feedback in Servo 
Systems* 


In addition to having a high degree of 
stability, a good servo system should have 
a quick response and no overshoot. But 
these two requirements are more or less con- 
tradictory. A simple analysis of a second or- 
der servo system shows that in order to have 
quick response, the frictional damping must 
be kept to a low value. But a low damping 
usually gives a large overshoot which can 
only be minimized by increasing the damp- 
ing. 

Lewis! in-his pioneering work showed 
that with a step input and a nonlinear feed- 
back proportional to the product of the 
“error” and the rate of change of output, it is 
possible to solve the problem. It was shown 
that such a nonlinear feedback results in an 
effective damping which is small at the be- 
ginning but becomes large as the output ap- 
proaches its final value. A low damping at 
the start gives a quick response. Again, no 
overshoot can occur, as the damping in- 
creases to a sufficiently large value before 
the output reaches its steady state level. 


The above type of feedback arrange- 


ment, which produces a very good result for 
a step input (for which it was analyzed), is 
not very effective if the input is of any other 
type. In the feedback arrangement of Lewis, 
the increase in damping is actually propor- 
tional to “error,” 7.e., difference between the 
input and the output. But as the input is 
constant for a “step,” the effective increase 
in damping becomes proportional to the 
output in this particular case. Thus, we get 
an increased damping as the output ap- 
proaches its steady state value and no over- 
shoot occurs. If, on the other hand, the in- 
put is not a constant step but itself changes 
with time as in a “ramp” or a “trapezoid” 
(step with finite slope), the effective increase 
in damping (as the output approaches its 
final value) will normally be small and over- 
shoots will occur. : 

If, however, we modify the system 
slightly by providing feedback signal pro- 
portional to the product of the output re- 
sponse and its time rate of change, better 


* Received by the PGAC, April 5, 1960. 
J. B. Lewis, Trans. AIEE, vol. 71, p. 449; 1953, 
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response is obtained for other types of in- 
puts as well. Such a feedback arrangement is 
discussed in the following paragraphs. 

In the servo arrangement discussed be- 
low, “G(s)” is the transfer function of the 
controller unit, and Ko is a constant giving 
amplifier gain; x; is the input signal to the 
system, and x is the resultant output. The 
first derivative (dx,/dt) of the output (xo) 
—as obtained from a rate generator—is 
multiplied with Koxo in a multiplier unit, 
and the product is added to the “error” 
e(=x;—4x.) in an adder. The resultant signal 
is then applied to the input of the controller 
unit and acts as the actuating signal for the 
controller. The arrangement is shown sche- 
matically in Fig. 1. Thus, the feedback sig- 
nal is composed of the sum of the error sig- 
nal and the product of output with its time 
rate of change. 


—— 


Fig. 1—Block diagram of the feedback arrangement. 


The differential equation of such a feed- 
back system may be represented by an equa- 
tion of the type 
dx 

at Dito +) = 0, 


ds i 
a2 ECR AS Bie) 


dt? 
where J and D are system constants, de- 
pending respectively on the inertia and the 
compliance of the system. R> is the frictional 
damping coefficient and B is a constant pro- 
portional to Ky. Or, we may write, 


x5 LR dXo 
dt* dt 


dj 


Jd + Dx = Dx; ++: (1) 


where 
R=R)+ Buo-::. (1a) 


Now, if we regard R as a constant over 
a small increment of time, for which the 
change in x, may be neglected, then the sys- 
tem can be treated as linear over that inter- 
val. Thus, the actual system is approxi- 
mated as a series of linear systems with in- 
creasing values of R. 

If R=R,(=Rot+Bxon) be the value of R 
during the mth increment of time, then dur- 
ing that particular interval 


Oxon dXon 
Rn 
di? i dt 


where %on is the value of x. at the mth inter- 
val. 

It is to be noted that (2) is identical with 
the equation derived by Lewis! if the input 
signal happens to be a step function (7.e., 
is equal to a constant). In that case, the 
above feedback arrangement will give iden- 
tical response with that of Lewis. The ef- 
fective value of R, in the Lewis system is 
given by 


if + Dion = Dxi(in) <--> (2) 


Rn. = Ro’ + Bgon — Xin) 


where x;,, is the value of x; at the mth instant. 
The above expression for R, becomes iden- 
tical with that of (1a) for the system under 
discussion only when x; is a constant, 7.e., 
for a step input. But if x; varies with time 
(as for all other types of inputs), the effec- 
tive damping in the above system will be dif- 
ferent from that of Lewis and is found to 
give a better response as shown later. 

Consider for example, a unit ramp input, 
for which x; =¢. In that case, (2) becomes 


FT AE i aa 
a AE Ey eg 


Solving this equation with the help of 
Laplace transform, we get 


aR (Ry l2I)t 
oe t —|— = 
c DLs D Re 
J 4]? 
R,2 D E, (ID) Rr? 
{Gi-2) (VFB) 
PRA fe J J 4]? 
Ry ID Ike 
Jakes PinboA J? 


o(V2 Ey. 0 


In (3), /D/J is the undamped natural fre- 
quency of the system and we may put 
V/D/J=. Also, for the case of critical 
damping, we have 


D oe 
(Rn)e = fap =2/JD--- (4) 


where (Rn) is the value of R, for the case 
of critical damping. 

It is evident from (1a) that the value of 
R, will increase with time as the output in- 
creases and the sinusoidal (and cosinusoidal) 
terms in (3) will die down rapidly with time, 
thus giving the steady state output of 
(t—R,/D) in a short while. But, due to low 
damping at the initial stage, the output will 
more closely follow the input than for the 
critically damped (or overdamped) case. A 
typical example is calculated and shown be- 
low. 

For the feedback arrangement under 
consideration, values of output response for 
the nth interval of time [as given by (3) | 
have been calculated successively by the 
method of numerical computation and suc- 
cessive approximation, as was also done by 
Lewis. In the example chosen to illustrate 
the result, the value of initial damping Ro 
was taken as 0.1(R,)- and the value of B was 
chosen as 0.9(R,).. Both the time (¢) and 
the output («) have been calculated and 
plotted in terms of wo. To calculate the value 
of R, we have taken wo as 6 radians per 
second. Any other value of wo would have 
given similar results. The calculated varia- 
tion of output with time for this typical 
case has been plotted in curve e of Fig. 2. 
In the same figure are also shown the char- 
acteristic curves for the input (curve a), the 
output with zero damping (curve 6), and 
the output with critical damping (curve c). 
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As expected, with zero damping, the output 
is oscillatory with large overshoots and 
with critical damping there is an appreciable 
time delay between the input signal and the 
output response. Both these results are un- 
satisfactory for most applications. But with 
the feedback arrangement shown schemati- 
cally in Fig. 1, there is no oscillation and no 
time delay (as shown by curve e). The over- 
shoot is also negligibly small. Steady state 
error is also found to be !ess than that under 
the critically damped condition. 

The output response with Lewis’s feed- 
back arrangement has also been plotted for 
comparison as curve d on the same figure 
(Fig. 2). Since, for the underdamped case, 
the error is never large and is actually zero 
at some points, there is no effective increase 
in damping with time. The output main- 
tains its damped oscillatory nature and 
appreciable overshoots occur. Thus, we see 
that with the feedback arrangement de- 
picted in Fig. 1, a step input gives the same 
type of improved response as that of Lewis, 
whereas for a ramp input a better response 
is obtained. 


Fig. 2—Output responses for the given system: curve 
b—output with zero damping; curve c—output 
with critical damping; curve d—output with Lewis’ 
feedback arrangement; curve e—output with the 
feedback arrangement under discussion (shown 
schematically in Fig. 1). The input ramp is shown 
as curve a for comparison. 


It may be mentioned that improved 
response occurs with the above feedback 
arrangement only if the initial damping of 
the system is low. But this is not a serious 
handicap. If the system is initially over- 
damped or critically damped, we can reduce 
the effective damping of the system by pro- 
viding an additional feedback loop with a 
suitable amount of positive feedback. 

B. CHATTERJEE 
Electrical Engrg. Dept. 
University of Illinois 
Urbana, III. 


Connection Between Various 
Methods of Investigating 
Absolute Stability* 

A general problem of control may be 
formulated as follows: Given a system (Fig. 


1) where f(c) may be any nonlinearity from 
the class defined by 


fle) = 0,0 = 0;. ofc) > 0, «0, (1) 


* Received by the PGAC, June 10, 1960. 
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what are the restrictions that must be placed 
on G to ensure absolute stability? Lure! 
discusses this problem in considerable de- 
tail and, using Lyapunof’s “second method,” 
derives after lengthy calculations rigorous 
stability criteria which guarantee absolute 
stability for any f(c) as defined above. 
These criteria and others derived by 
Lyapunef’s “second method”? are difficult 
to apply to complicated systems, and phys- 
ical insight into the underlying reasons for 
stability is often lost because of the cumber- 
some calculations. 

Another approach which clarifies this 
problem is given by Sobolev,? who compares 
the nonlinear system (Fig. 1) to a linear 
system (Fig. 2). He proves the following 
theorem: 

If the linear system is stable for all 
0<h<o«, then the nonlinear system is 
absolutely stable. Thus, absolute stability of 
the nonlinear system is given by the require- 
ment for unconditional stability of the as- 
sociated linear system. Hence, the equation 


1+4G=0, 0<h<o (2) 


must have no roots with positive real parts. 


Fig. 2. 


Sobolev solves this equation by using 
the familiar Routh-Hurwitz criteria. The re- 
sulting equations define the stability region 
of the linear system with # as parameter. 
The envelope of this function gives the sta- 
bility region of the linear system for all h, 
and hence by the above theorem the region 
of absolute stability of the nonlinear sys- 
tem. 

The above stability condition (2) can 
also be written as 


Oe ieols (3) 


Hence, the stability region for the nonlinear 
system is given by 


/G(jw) > — 180°. (4) 


In this form, the significance of Sobolev’s 
theorem is made clearer. Not only because 
the calculations are simpler, obviating the 
need for finding the envelope, but also, 


1A. 1, Lure, “Some Non-Linear Problems in the 
je, of Control,” H. M. Stationary Office, London; 
195 


2 A, K. Bedel’baev, “Some simplified stability cri- 
teria for non-linear controlled systems.” Azxtomat. 
i Telemekh., vol. 20, pp. 689-701; June, 1959. 

3 Yu. S. Sobolev, “The absolute stabilities of cer- 
tain controlled systems,” Automat. 1 Telemekh., vol. 
20, pp. 401-405; April, 1959. 
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mainly, because the role of the transfer func- 
tion in determining the stability of the non- 
linear system is placed very much in evi- 
dence. 

Eq. (3) has a striking resemblance to the 
equation for determining stability by the 
describing function method. 


1 

. N(o) ©) 
where JV(c) is the describing function for the 
nonlinearity f(c). By the above assumptions 
on f(¢) (1), N(c) is real and positive. Hence, 
by describing function analysis, stability in 
the large will again be ensured only if 
/G(jw) > —180° for all w. Thus, the approxi- 
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mate describing function method yields the ae 
same stability criterion as given by Sobolev’s 
rigorous approach. - a>o 
Fig. (4) also enables a simpler approach b 
to the calculation of the boundaries of the , 
stability region. Consider Sobolev’s first 4 Sl de 
example (Fig. 3). de 
The transfer function is J 
Vs 
a 1 [ i Nes Es+ —] 
s(W + T;?s) s(U + T,?s) 
= 1 s?(NG?+ Ta*r) + s(ur + NE) + Na Fac 
2 2T2T,? + s(uT 2 + WT?) + uW Fig. 5 
1 s*M,+5M;+ Na A d \ pe Re 
=—: = s a second example, consider Lure’s 
3? °T?T.? + sMi + UW general canonical form, 
1 P(s) Lp = Lp tf (e), p=1,2,---,n 
=—: ‘ ‘ 
s? Q(s) «=> BZ) — rf(o). (7) 
The stability boundary is given by = 
By eliminating all Z,’s, we obtain the open 
/G(jw) = — 180°. (6) loop transfer function 
ae ie ay 
=—\|r ——— |. 
Hence, S ae geen 
/P(je) — /Q(jo) “ee For n=2, the transfer function becomes 
G 1 5% — s[r(\r + 2) + Bi + Bo] + rrid2 + Bid2 + B21 
ae tg 5? — s(Xi + Az) + Arde 
or applying again the stability boundary cri- 
terion G(jw) = 180°, we arrive after some al- 
ia oMs 2 oMi a gebra at the bi-quadratic equation 
tan—1 ———__—__ tan 1 ———__-___ = 0 
Na—w*M2 UW — wT ?T 2” 


After some algebra, we arrive at 
/ uW M3 faa Na: My, : 
°" V T2T 2M; — MoM, 


The system will be stable only if this 
equation has no real roots for all w. The 
stability region is shown in Fig. 4. Substi- 
tuting for Mi, Ms, and M;, we arrive at 
Sobolev’s conditions: 


UW ur + NE) — Na(uT;? + WT.”) = 0 


NT#(ET2? — Gu) — (NG? + Ta*r)WTd? < 0. 


ctr + wo?[r (da? + Ao?) + Bid + Bare] 


Bi , Be 


+)1?do? [r++ | =w!-a+w?-b +6 
1 


a5 
= 0, 
The stability region in this case is shown 
in Fig. 5. 
One boundary is given by c20, or 

Bi , Be 
= — + — > 0(M1, Ax ¥ 0). 8 
T WE Se ign (A1, Az ) (8) 


This is Lur’e’s first condition.t The sec- 
ond boundary is given by the lower branch 
of the parabola 6?=4ac, or 


4 Lure, op. cit., p. 71, (9.6). 
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[(Ar? + do?) + Bids + Bode]? 
= Ard2r2T?, (9) 
Lure’s second condition is* 
mr? = rf? — vf, + By + DevV/rTf2 = 0. (10) 
Upon substituting 
e= 11; By = Biri + Baro; 
Fi = M1 + Yo; fo = Arde 


We see that conditions (9) and (10) are 
the same. Hence, the simple criterion 
G(jw) > —180° leads to the same results as 
obtained by Lyapunof’s second method. 

In conclusion,® Sobolev’s comparison the- 
orem leads to a simple and rigorous stability 
criterion. Also, it gives a firm basis to stabil- 
ity analysis by the describing function 
method. 

ZE’EV BONENN 

Dept. of Engrg. 
Cambridge University 
Cambridge, England 

5 Note added in proof: After writing this corres- 
pondence the author was told during the first IFAC 
Congress in Moscow that Sobolev’s theorem is not 
always valid for n >2. A lively correspondence on this 
subject has appeared in recent issues of Aviomatika 1 
Telemekhanika, It has not yet been possible to obtain 
an English translation of these issues. However, the 


results of the above examples are correct and agree 
with those derived by Lyanimof’s second method. 


On the Approximation of Roots of 
nth Order Polynomials* 


The roots of an nth order polynomial 
may be obtained by a wide variety of tech- 
niques. Lin’s method! of quadratic factor 
extraction is particularly useful, because of 
its simplicity, when one is dealing with poly- 
nomials having real coefficients where it is 
known that complex roots will occur in con- 
jugate pairs. While Lin’s method converges 
to one of the quadratic factors of the poly- 
nomial, convergence is sometimes slow and, 
therefore, hand computation of a reasonable 
approximation can be tedious. 

The root-locus method? provides an ap- 
proximation which may then be used as a 
first trial factor in Lin’s method for final 
convergence. This technique requires that 
the polynomial be rearranged into the fa- 
miliar root-locus form, 


KN(s) _ 
Dis) s 


FS (1) 


One convenient way of accomplishing the 
rearrangement is to divide the polynomial by 
all but the three lowest order terms as indi- 
cated as follows: 


* Received by the PGAC, June 1, 1960. The devel- 
opment discussed here was accomplished while the 
author was employed by the Microwave Res. Inst., 
Polytechnic Inst. of Brooklyn, Brooklyn, N. Y : 

1S. N. Lin, “A method of successive approxima- 
tions of evaluating the real and complex roots of cubic 
and higher order equations,” J. Math. Phys., vol. 
20, pp. 231-242; August, 1941. : 

2 J, G. Truxal, “Root locus methods,” in “Auto- 
matic Feedback Control System Synthesis,” pp. 267— 
272, McGraw-Hill Book Co., Inc., New York, N. Y.; 
1955, 
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P(S) = 5” + ays"? + das” + + + © + Gy28" + dn15 + On = 0 (2) ie 
An—-1 Qn 
Gn—2 (3 + —s+ =) 
Gn—2 An—2 

=1 2 : 
0) os Se(sra8 bE aysn—4 + dos75 + Bar. ate pe ( ) 

2(s? nie fl 
O=1+4 An—2(s Sm 15 => Gy) K 


If d,_2 is now considered a variable, (4) 
is of the same form as (1). With x less than 
six, factors of NV(s) and D(s) will, at most, be 
quadratics and their roots may be easily 
determined, thus locating pole-zero posi- 
tions. In such a case, the locus of pole posi- 
tions as K varies may be plotted directly. 
The locations of poles for K =dn_2 are the n 
roots of the original polynomial in s. With x 
six Or more, a repeated application of this 
technique is suggested. That is to say, 


D(s) = s3(s"-8 + ays™“4 + +++ +n) (5) 


is set equal to zero and its roots are deter- 
mined by means of the root-locus technique. 

A useful extension of the root-locus tech- 
nique is obtained if one notes that more than 
one representation for the root locus can be 
obtained from the same polynomial. If both 
loci are plotted, they must intersect at the 
roots of the original polynomial. Consider, 
for example, the following fourth order poly- 
nomial: 


s¢+as?+ bs? +estd=0. (6) 


Dividing by all but the three lowest order 


terms yields 
b G ae = Se +) 
b b 
s3(s + a) 


Division by all but the two lowest order 


terms yields 
d 
(s+7) 
C 


Cennard ge B) (8) 


Uh thes (7) 


The intersecting root-loci technique de- 
scribed in these notes may be applied to 
polynomials of order 5 or greater by repeated 
application of the method to reduce denom- 
inator factors to quadratics. 

It should be noted that loci will intersect 
over a range of values for real roots and, in 
general, at a point only for complex roots. 
Given such a range of values for real roots, 
it is often possible to reduce the order of the 
polynomial by trial-and-error removal of the 
real roots or by means of Lin’s method ap- 
plied to single real roots. 

By way of illustration, consider the 
fourth-order polynomial 


P(s) = st 15s + 55s? + 245s + 204=0. (9) 


Arranging the polynomial in a root-locus 
form, 


245 204 
ol Cara 
ee ST Re 


"which reduces to 


K,i(s + 1.09) (s + 3.41) 


s3(s + 15) pare ee 


es 


SOS Oe ars © fe gs 8 So Eg, 3) 


A second root locus form is obtained as 


follows: 
245 ( of ss) 
te 48 =0 (12) 
sX(s?@+155+55) ’ 
K2(s + 0.83) 
P35 52(s + 6.4)(s + 8.6) _ o. 2) 


The loci for variable K in (11) and (13) 
are sketched in Figs. 1 and 2. Fig. 3 is a 
sketch of the intersections of the two loci, 
and thus indicates the roots of the poly- 
nomial of (9). 

Examination of Fig. 3 indicates a real 
root in the range —0.83>s>—1.09. One 
might well assume the root to be located at 
s=—1 (as, indeed, it is) and thereby reduce 
the polynomial. The reduced polynomial, 


ERs age aiy 04 Oa) 
(s + 1) ; 
may be written in three root-loci forms: 
K 4.98 
ee, (15) 
s?(s + 14) 


: Kale + 1.87 +j3.48) (s+ 1.87— 3.48) 


; (16) 


53 
and 
te us 
s(s + 4.14)(s + 9.86) 


A sketch the intersections of these loci is 
presented in Fig. 4. 

On comparison of the intersecting root- 
locus technique with Lin’s method of quad- 
ratic factor extraction, it is noted that, for 
the polynomial of (9), an oscillatory solution 
results from Lin’s method, while Fig. 3 


0. (17) 


Fig. 1—Root locus sketch for 
Ki(s + 1.09)(s + 3.41) 


s3(s + 15) 


Fig. 2—Root locus sketch for 
K2(s + 0.83) 


+ = 0. 
s*(s + 6.4)(s + 8.6) 


fj 


| 
| 


Fig. 3—Intersecting loci sketch for 
P(s) = st + 1553 + 555? + 245s + 204 = 0 
= (s +1)(s + 12)(s? + 2s +17) =0. 


Fig. 4—Intersecting loci sketch for 
P(s) = 5s? + 14s? + 415 + 204 =0 
= (s +12)(s2 + 2s +17) =0. 


1960 


yields a close approximation to the complex 
roots and one of the real roots. If one applies 
Lin’s method to remove a real root from 
(9), convergence to the factor (s+1) is 
obtained in three trials and the polynomial is 
reduced by an order. Reapplying Lin’s meth- 
od to obtain the complex roots (or the re- 
maining real root), the complex roots are 
approximated by the third or fourth trial 
divisor. 

The described extension enhances the 
value of the basic technique by indicating 
the regions of the roots from crude sketches 
of the root-loci. Increased accuracy may be 
obtained by any of three methods: 


1) Both loci may be plotted carefully in 
the regions of intersection. 

2) One of the loci may be plotted ac- 
curately in the regions of intersections 
and the points corresponding to K 
located on the locus. 

3) The approximate intersection ob- 
tained from the crude sketch may be 
used as a first trial factor in another 
technique such as Lin’s or Newton's 
method. 

Joun D. GLomMB 

Central Res. Engrg. Div. 
Continental Can Co. 
Chicago, Ill. 


Determination of Transfer Func- 
tion Coefficients of a Linear 
Dynamic System from Frequency 
Response Characteristics* 


Determination of transfer functions in 
standard form, from experimentally ob- 
tained frequency response characteristics, is 
a classical problem of extraordinary interest 
to designers of servo-systems. Levy! has re- 
cently presented a very versatile method 
which facilitates exact synthesis of frequency 
dependent systems. The method can be used 
to fix up an algebraic expression as the ratio 
of two polynomials, which could later be 
solved, with the aid of an Isograph,? to give 
the zeroes and poles of the system under 
study. 

However, as pointed out by Levy him- 
self, the utility of his method is limited by 


two major restrictions. The object of this, 


correspondence is to suggest means of over- 
coming these restrictions in a very direct 
and elegant manner. The first restriction ac- 
cording to Levy is that systems with infinite 
gain at zero frequency cannot be handled by 
his method except by an indirect approach. 
These sytems are the familiar type I, type 
II, and type III Servo-systems. It is inter- 
esting to observe that this restriction is 
totally removed if one uses M-curves instead 
of G-curves. G and M refer to the open and 
closed-loop responses in the frequency do- 
main, as per standard notation. M-curves 


# Received by the PGAC, June 8, 1960. 
1 ag ale pad “Complex curve fitting,” IRE 
TRANS. ON Automatic ContRoL, vol. AC-4, pp. 
—43; May, 1959, : 
oH 7B. Weakata Rao, “A novel type of isograph,” 
IRE TRANS. ON ELECTRONIC COMPUTERS, vol. EC-7, 
pp. 97-103; June, 1958. 


tit ie : a 


Correspondence 


can be easily obtained from the G-curves 
with the aid of a Nichols Chart. The use of 
M-curves, in preference to G-curves, results 
in several additional advantages. Firstly, in 
the case of unity feedback systems the de- 
nominator polynomial of the M-curve di- 
rectly gives the characteristic equation of 
the system, which is useful for determining 
quantitatively the performance character- 
istic of the system in the time domain. Sec- 
ondly, computation with M-curves is much 
easier because of the smooth nature of these 
curves. In any event, if one wishes to get the 
transfer function for G, it can be obtained as 
follows. 

Let 


fiGjo) 
fa( jo) 


where jf; and fz are polynomials in jw. As per 
standard notation, 


M(je) = 


(1) 


AGES 
M (jw) = 1+ GGo) (2) 

From (1) and (2), 
oa 6) 


foljo) — faljoo) — 


The necessity for using M-curves in case 
of type I, type II or type III systems can be 
illustrated by a specific example. 


Let 
G(jo) =f (jw): ; » for n= 1; 2,3, ae) (4) 
(jw)” 
where f(jw) is of type 
Ao + Ai(jo) + Aa(Jo)? +o (5) 


1 + Bi(jw) + Bo(jw)? + +-- 


Eq. (4) obviously gives infinite gain at zero 
frequency. Considering the corresponding 
M-function, one gets from (2) and (4), 


epee) 
Goo)” + fa) 
Eq. (6) gives a finite gain even at zero fre- 
quency, thus overcoming the first restriction. 
The second restriction of Levy’s method 
pertains to the weighting of the function for 
the error criterion. The use of the type of 
weighted error function suggested by Levy 
implies that for a given value of w, the mag- 
nitude of the error is directly proportional 
to the magnitude of the function. In order 
to minimize the errors in computation, Levy 
sugggests selecting a greater number of 
sampling points in the critical region of the 
curve. A more direct means of overcoming 
this restriction is by working in terms of 
the inverse plots. By working with the in- 
verse function, the weighting is such as to 
make the error a minimum in the neighbor- 
hood of the local maxima of the given func- 
tion. 
eet 


M (jw) (6) 


, (7) 


be the required function, where 7(jw) is the 
numerator and d(jw) the denominator of 
M(jw). Let 


H(jo) = Rw) + JI@) (8) 


be the experimental curve, with R(w) denot- 


ing its real part, and J(w) its imaginary part. 
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Let the numerical difference between the 
two functions 1/M(jw) and 1/H(jw) repre- 
sent the error in curve-fitting, that is 

il d() 
H(jw) nw) 
Multiplying both sides of (9) by n(w), one 
gets 


5(@) = (9) 


n(w) 
H (jw) 
The RHS of (10) is a function of real and 


imaginary terms, which may be separated to 
give 


5(@)-n(w) = — d(w). (10) 


5(@) nw) = aw) +76); (11) 


or at any specific value of frequency w=wx, 


[5(w)-m(@) Joox = a (wg) +8 (wx). (12) 


Now, the weighted error E is defined as the 
function given in (9) summed over the sam- 
pling freqencies wx. Thus, 
E= Di [ar(og) + 0%(eg)]. (13) 
The unknown coefficients are determined 
on the basis of minimizing the function F. 
By the above derivation, it is clear that 
[5(w) -n(w)]? is a minimum, say K2, for that 
set of coefficients determined by the princi- 
ple of least squares of the weighted eiror 
function. For this choice, 6(w) is a minimum 
K/n(), where n(), or the numerator of 
M(jw) is a maximum. Thus, the error in 
fitting is a minimum in the neighborhood of 
local maxima of M, an aspect of obvious 
advantage for design purposes. It is easily 
seen that the suggested modifications widen 
the scope of Levy’s method to a considerable 
extent without introducing any new errors. 
Further, the constant term of the nu- 
merator polynomial is equal to the magni- 
tude of the function at negligibly small val- 


‘ues of w. Utilizing this data reduces the 


order of the matrix equation by one. 

Recently, Dudnikov? has suggested an 
interesting method of determining the co- 
efficients of the polynomials from the initial 
portion of an experimentally obtained am- 
plitude-phase characteristic. 

The transformer analog principle* can be 
advantageously applied for developing a 
special purpose computer to mechanize the 
the evaluation of the manipulated parame- 
ters \, S, Zand U. It is proposed to use this 
equipment in conjunction with a Mallock’s 
equation solver,® to determine straightaway 
the coefficients of the numerator and de- 
nominator polynomials. Work on the de- 
velopment and design of this special purpose 
analog computer is under way. 

The authors wish to express their sincere 
thanks to Prof. H. N. Ramachandra Rao 
for his kind interest in this work. 

P. VENKATA Rao 
(Mrs.) C. LAKsumi Bat 
Dept. of Power Engrg. 

Indian Inst. of Sci. 
Bangalore, India 


3E. E, Dudnikoy, “Determination of transfer 
function coefficients of a linear system,” Automation 
and Remote Control (English translation of Avtomat. ¢ 
Telemekh.), vol. 20, pp. 552-559; May, 1959. 

4P, Venkata Rao and G. Krishna, “The trans- 
former analogue computer (TAC),” Trans. AIEE, 
vol. 34, pp. 732-739; January, 1958. 

5M. W. Humphrey Davies and G. R. Slemon, 
“Transformer analogue network analyzer,” Proc, IEE, 
vol. 100, pt. 2, pp. 469-486; 1953, 
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Notice 
INCREASE IN GROUP MEMBERSHIP FEE 


The PGAC is endeavoring to provide the greatest 
possible service to its members in the publishing of 
control papers. The rate of submission of good TRANS- 
ACTIONS papers continues to increase and we are nearing 
our goal of a quarterly TRANSACTIONS on a regular 
schedule. In addition, the PGAC plans to continue the 
practice started in 1959 of publishing the PGAC por- 
tions of the IRE INTERNATIONAL CONVENTION and 
WESCON REcorDs as special TRANSACTIONS issues, 
and will also publish IRE papers presented at the an- 
nual Joint Automatic Control Conferences. Needless to 
say, publication of the PGAC TRANSACTIONS costs 
money, and the PGAC must keep a balanced budget. 

Each of the Professional Groups has its own funds, 
derived in part from an IRE subsidy. For a number 
of years the IRE provided support to the Professional 
Groups by matching membership fees, but this was 
recently changed in a manner to emphasize and promote 
Group publications. Thus, in 1959 the IRE supported 
the Professional Groups by providing one-third of all 
costs directly attributable to the publication of TRANS- 
ACTIONS and Newsletters. However, the rapid increase 
in publications (and therefore in IRE subsidies) planned 
by the 28 Professional Groups has forced the IRE Board 


of Directors to depart from strict application of the one 
third rule in 1960, and to place ceilings on the support 
provided each group. Thus, the IRE publications sub- 
sidy to the PGAC in 1960 will be limited to $3000, 
rather than the figure of $4000 submitted in the PGAC 
budget. 

With the reduced level of IRE support, the PGAC 
balance sheet shows a negative derivative which ob- 
viously cannot continue for very long. Since it is not 
anticipated that the level of [RE support will be signifi- 
cantly increased in subsequent years, some action was 
clearly required. The PGAC Administrative Com- 
mittee discussed this situation at its meeting in March, 
1960, and decided to raise the PGAC membership fee 
from $2.00 to $3.00, effective July 1, 1960. The only 
alternative, the reduction of PGAC publications, was 
not considered desirable. It is noted that of the 28 IRE 
Professional Groups, 11 already have a $3.00 fee, and 2 
have a $4.00 fee. 

It is hoped that this decision of the PGAC Adminis- 
trative Committee to continue present publication 
policies will meet with the approval of the PGAC mem- 
bership. 

Joun E. Warp, Chairman, PGAC 
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Hazeltine, Barrett 
Hersh, A. I. 
Heuchling, T. P. 
ie BO De oo oe 

Hills, F. B, 

Hills, W. L. 

Ho, Yu-Chi 
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Johnson, E. P., Jr. 
Kailath, Thomas 


Keith, G. E, 
Kelleher, J. J. 
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Olsson, E, A. 
Orenberg, Arthur 
Ortolano, A. J. 
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Parke, N. G., III 
Parkinson, B. N. 
Passera, A. L. 
Pastan, H. L. 
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Rousseau, W. F. 
Rowell, W. G. 
Roy, A. G., Jr. 
Sabin, E. A. 
Saunders, N. B. 
Scanlon, W. C. 
Schnall, Emanuel 
Schochet, J. R. 


Schramm, M, W., Jr. 


Schroeder, R. L. 
Schwerdt, K. R. 
Seifert, W. W. 
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Shappy, R. A. 
Sheingold, D. H. 
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Wolf, S. H 


Wolfe, Russell 
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Jones, D. R. 
Lightburn, R. C, 
Martin, F. C. 
Meister, L. E. 
Merle, C. W. 
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Powell, F. D. 
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Rynaski, E. G. 
Schneeberger, R. F. 
Schultz, W. C. 
Stieber, Alexander 
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Bourret, C. J. 
Brockett, R. I. 
Brooks, F. A., Jr. 
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Mapes, T. J. 
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Dodson, G. C., Jr. 
Ellert, F. J. 
Kirchmayer, L. K, 
Lippitt, D. L. 
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1960 


Ashcroft, D. L. 
Bessette, D. U. 
Brady, D. J. 
Brule, J. D. 
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Burke, T. H. 
Clark, P. B. 
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Croly, J. W. 
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Schulkind, Donald 
Schuster, G. M. 
Schwartz, Albert 
Schwartz, Bernard 
Schwartz, I. S. 
Scott, J. E. 
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Winzemer, A. M. 
Witmer, F. S. 
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